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Abstract 

The  United  States  Air  Force  is  interested  in  the  potential  side  effects — at  the 
cellular  level — from  exposure  to  mission-essential  chemicals.  Presently,  Air  Force 
toxicology  studies  are  conducted  to  help  shed  light  in  identifying  potential  hazards 
to  workers.  However,  it  takes  a  considerable  amount  of  money,  resources,  and  time 
to  obtain  and  analyze  experimental  results  from  toxicology  studies.  The  necessity  for 
innovative  methods  that  enable  researchers  to  more  effectively  generate  and  analyze 
data  is  apparent. 

Mathematical  modeling  is  a  viable  option  to  become  a  valuable  tool  for  the 
researcher.  Mathematical  models  can  rapidly  generate  informative  predictions  on 
how  a  cell  reacts  to  a  certain  toxicant  exposure.  Moreover,  information  is  readily 
available  when  generated  by  mathematical  models. 

This  research  involves  the  study  of  one  non-biological  reaction  system  and  four 
biological,  intracellular  reaction  systems.  Each  system  is  converted  into  a  mathe¬ 
matical  model  using  the  rate-equation  approach.  Numerical  simulation  results  from 
these  mathematical  models  are  obtained  using  two  novel  software  modeling  tools 
and  MATLAB.  Results  obtained  from  the  novel  modeling  tools  are  compared  to 
MATLAB’s  results  in  order  to  ascertain  the  accuracy  of  each  novel  modeling  tool. 

The  experience  that  is  gained  in  deriving  mathematical  models  and  using  novel 
tools  to  perform  numerical  simulations  for  these  reaction  systems  should  help  the 
Air  Force  develop  intracellular  models  to  assist  in  future  toxicology  studies. 


DETERMINISTIC  INTRACELLULAR  MODELING 


I.  Introduction 

1.1  Overview 

The  potential  of  Air  Force  personnel  being  exposed  to  harmful  chemicals  within 
the  work  environment  must  be  taken  under  serious  consideration.  This  risk  to  per¬ 
sonnel  is  heightened,  because  Air  Force  members  are  oftentimes  required  to  work — 
over  a  wide  spectrum  of  conditions — on  or  around  specialized  equipment  that  contain 
an  assortment  of  toxic  substances.  Toxicity  studies  are  conducted  to  help  shed  light 
in  identifying  potential  hazards  to  workers.  However,  the  cost  of  toxicity  studies  is 
high. 

Performing  toxicity  studies  are  formidable  for  two  reasons.  Experiments  are 
very  costly  in  funds  and  resources.  Also,  it  usually  takes  a  considerable  amount 
of  time  to  obtain  and  analyze  experimental  results.  The  necessity  for  innovative 
methods  that  enable  researchers  to  more  effectively  generate  and  analyze  data  is 
apparent. 

Mathematical  modeling  has  the  potential  to  become  a  valuable  tool  for  the 
researcher.  Credible  mathematical  models  can  generate  informative  predictions  on 
how  a  cell  reacts  to  being  exposed  to  a  certain  toxicant.  Moreover,  information  is 
readily  available  when  produced  by  mathematical  models. 

1.2  Problem 

In  order  for  the  Air  Force  to  develop  mathematical  models  that  can  be  ulti¬ 
mately  coupled  with  existing  experimental  data  on  the  cell,  it  is  important  to  under¬ 
stand  how  various  intracellular  processes  are  regulated.  By  deriving  computational 
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models  from  well-documented  cellular  systems  and  implementing  those  models  on 
novel  software  tools,  valuable  experience  is  gained.  This  newfound  knowledge  can 
be  used  to  model  cellular  systems  of  a  higher  order.  This  knowledge  will  lead  to  a 
better  understanding  of  how  the  Air  Force  should  construct  intracelluar  models  to 
assist  future  toxicity  studies. 

1.3  Scope 

This  research  involves  a  literature  review  on  biological  topics  relevant  to  the 
living  cell.  This  review  includes  microbiology,  biochemistry,  chemical  kinetics,  and 
intracelluar  structure.  Mathematical  models  (deterministic)  are  developed  from  syn¬ 
thetic  and  naturally  occurring  networks.  These  models  are  subsequently  simulated 
on  three  software  applications:  MATLAB,  BioCharon  ([2]  and  [5]),  and  JigCcll  ([31] 
and  [30]).  Results  from  these  software  applications  are  analyzed. 

1-4  Summary  of  Thesis 

This  thesis  is  organized  as  follows: 

Chapter  2  gives  an  overview  on  cellular  biology  and  biochemistry.  Then,  math¬ 
ematical  modeling  is  formally  introduced.  Three  alternative  approaches  for  mathe¬ 
matical  modeling  are  briefly  discussed.  It  concludes  with  a  discussion  on  exisiting 
software  tools  that  can  be  used  for  intracelluar  modeling. 

Chapter  3  gives  a  detailed  description  of  the  rate-equation  approach  (a  math¬ 
ematical  modeling  method).  Three  software  applications — MATLAB,  JigCcll,  and 
BioCharon — that  can  be  used  to  model  intracelluar  processes  are  briefly  discussed. 
One  non-biological  network  and  four  naturally  occurring  networks  are  presented. 
Computational  models  for  the  synthetic  and  three  naturally  occurring  networks  are 
explicity  derived. 

Chapter  4  presents  general  simulation  results  and  comparisons  for  the  reaction 
systems  discussed  in  chapter  3.  Both  general  and  comparison  results  are  analyzed. 


1-2 


Chapter  5  summarizes  the  work  completed.  It  gives  conclusions  and  recom¬ 


mendations  for  future  work. 
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II.  Background 


2. 1  Overview 

“Cells  are  the  fundamental  working  units  of  every  living  system  [29:1].”  They 
encase  a  dynamic  environment  of  ongoing  activities,  such  as  synthesis  of  proteins 
or  the  breakdown  of  glucose.  The  underlying  logic  for  cellular  function  is  highly 
complex.  Understanding  this  functionality  is  the  bedrock  for  conceptualizing  and 
comprehending  higher  order  phenomena,  such  as  physiology,  anatomy,  and  ecol¬ 
ogy  [14:544], 

The  dynamic  environment  within  the  cell  (intracelluar)  involves  a  sophisti¬ 
cated  degree  of  interaction  between  three  important  classes  of  macromolecules:  de¬ 
oxyribonucleic  acid  (DNA),  ribonucleic  acid  (RNA),  and  proteins  [14:544],  These 
molecules  continuously  interact  by  means  of  biochemical  reactions.  Some  of  these  re¬ 
actions  take  place  independently  and  others  take  place  serially.  As  a  result,  molecules 
within  the  cell  are  in  constant  flux.  This  flux  is  oftentimes  regulated  by  genes. 

Genes — specified  segments  of  DNA — control  intracellular  processes  by  spec¬ 
ifying  the  synthesis  of  enzymes  (a  special  form  of  a  protein)  and  other  types  of 
proteins  [9:6].  Specific  groups  of  genes  may  be  activated  or  ‘expressed’  by  particular 
signals  in  order  to  regulate  a  common  cell  process.  Moreover,  groups  of  genes  may 
modulate  expression  levels  of  other  genes.  Such  groups  are  called  genetic  regulatory 
networks  [26:248]. 

These  types  of  networks  are  commonly  divided  into  separate  modules  of  opera¬ 
tions.  In  this  way,  cell  functionality  can  be  viewed  “as  a  collection  of  interrelated  sub¬ 
systems  [14:544].”  Nonetheless,  describing  these  networks — especially  quantitatively — 
is  a  monumental  undertaking.  Smolen  et  al.  elaborate  on  the  inherent  complexity 
of  genetic  regulatory  networks:  “Understanding  the  combined  effects  of  these  phe¬ 
nomena  is  often  beyond  the  capacity  of  intuition  [26:249].” 
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A  common  approach  to  begin  the  elucidation  of  cellular  functionality  is  through 
modeling  methods.  Models  describing  cellular  networks  must  embody  hypotheses 
made  about  each  subsystem  as  well  as  on  how  these  subsystems  are  subsequently 
integrated  back  together.  Presently,  mathematical  modeling — via  mapping  a  cellular 
network  to  a  system  of  ordinary  differential  equations  (ODEs) — is  widely  used  to 
describe  the  dynamics  of  a  cellular  network.  Hasty  et  al.  describe  the  process  of 
implementing  mathematical  models  [10:277]: 

The  modeling  of  gene  regulatory  networks  relies  on  characterization  of  the 
behavior  of  small  subsystems,  formation  of  hypotheses  about  how  these 
subsystems  interconnect,  translation  of  these  hypotheses  into  a  mathe¬ 
matical  model  and  experimention  to  yield  results  that  indicate  necessary 
changes  to  the  original  hypotheses. 

Mathematical  modeling  has  two  major  advantages.  “The  precision  of  the  math¬ 
ematical  language  makes  mathematical  modeling  a  useful  framework  for  conceptual¬ 
izing  and  understanding  complex  biochemical  systems  [26:248].”  Consequently,  the 
modeler  is  able  to  precisely  define  hypotheses  made  about  each  network  subsystem 
and  overall  network  organization  by  a  well-defined  language.  Most  importantly, 
present  circumstances  (explained  below)  are  primed  for  quantitative  as  well  as  qual¬ 
itative  coupling  to  experimental  data. 

Since  large  amounts  of  data  at  the  genomic  level  are  readily  available  and  exper¬ 
imental  methods  are  continually  being  enhanced,  the  stage  is  set  for  building  models 
that  make  physical  sense.  Consequently,  the  modeler  is  able  to  ‘verify’  candidate 
models.  Plus,  the  modeler  can  make  ‘informed’  modifications  to  model  parameters 
(e.g.,  change  values  of  kinetic  rate  constants),  based  on  experimental  data.  In  this 
way,  mathematical  models  have  the  potential  to  be  tightly  coupled  to  experimental 
data. 

For  example,  the  modeler  can  compare  a  temporal  data  set  for  a  natural  bio¬ 
logical  system  to  model  predictions.  Discrepancies  between  the  data  set  and  model 
predictions  are  a  strong  indication  that  the  proposed  model  contains  errors.  The 
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modeler  can  then  refine  or  even  reformulate  certain  hypotheses  made  about  net¬ 
work  architecture  and  exactly  implement  those  changes  in  the  mathematical  model. 
Since  formulating  a  computational  model  involves  a  series  of  intricate  steps — often 
overwhelming  the  pencil-and-paper  method — there  is  a  necessity  to  develop  software 
packages  that  assist  modelers  in  the  modeling  process. 

Fortunately,  there  are  already  existing  software  packages  for  modeling  biochem¬ 
ical  networks.  In  general,  these  software  packages  add  a  beneficial  layer  of  abstraction 
between  a  schematic  or  wiring  diagram  that  describes  the  architecture  or  pathways 
of  a  biochemical  network  and  the  the  corresponding  computational  model  generated 
from  that  diagram.  This  layer  of  abstraction  is  often  incorporated  by  means  of  a 
graphical  user  interface  (GUI).  The  user  enters  details  about  the  biochemical  net¬ 
work  in  symbolic  form  (e.g.,  graphical  objects  or  chemical  reaction  notation).  Then, 
the  software  generates  the  computational  model  from  this  input.  In  this  fashion,  the 
user  is  spared  the  intracacies  of  explicitly  transforming  (mapping)  the  schematic  or 
wiring  diagram  to  its  corresponding  computational  model.  Consquently,  the  likeli¬ 
hood  of  making  errors  in  deriving  compuational  models  is  usually  reduced,  especially 
when  attempting  to  model  large  networks. 

As  a  precursor  for  using  available  software  packages,  the  modeler  must  have  an 
understanding  of  the  fundamentals  of  life  sciences  related  to  cellular  networks,  such 
as  microbiology,  biochemistry,  and  genetics.  Since  mathematical  modeling  generally 
requires  derivation  of  a  system  of  nonlinear  ordinary  differential  equations  (ODEs) 
from  a  wiring  diagram,  the  modeler  should  also  be  well-grounded  in  mathematical 
analysis  of  nonlinear  ODEs.  From  a  mathematical  (technical)  perspective,  the  mod¬ 
eler  is  now  able  to  predict  and  explain  subsequent  model  behavior.  In  the  following 
sections,  I  briefly  elaborate  on  these  essential  topics  that  serve  as  a  prerequisite  for 
modeling  cellular  networks. 
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2.2  Cellular  Biology 

Understanding  cellular  structure  is  essential  in  implementing  the  correct  ap¬ 
proach  to  model  cellular  behavior.  Cellular  functionality  is  generally  partioned  into 
various  components  of  the  cell.  In  order  to  model  mechanisms  of  cellular  behavior, 
these  cellular  components  must  be  understood.  This  section  contains  references  to 
intracellular  structure  for  prokaryotic  and  eukaryotic  cell  types. 

2.2.1  Cell  Types.  All  living  cell  types  can  be  classified  as  either  prokaryotic 
or  eukaryotic  [6:71].  The  structural  complexity  of  prokaryotes  is  simpler  than  eu¬ 
karyotes.  Organisms  in  the  prokaryotic  class  consist  of  only  one  cell,  such  as  bacteria. 
They  “lack  a  nucleus  and  other  membrane-enclosed  structures  [6:71].”  On  the  other 
hand,  eukaryotes  encompass  all  plants,  animal,  fungi  and  protists  [6:71].  Structures 
in  this  class  are  more  defined.  For  example,  cells  in  this  class  possess  a  nucleus 
and  other  membrane-enclosed  structures.  Related  material  about  prokaryotic  and 
eukaryotic  cell  types  can  be  found  in  [6:71-95]. 

2.2.2  Intracellular  Structure.  Using  bacteria  cells  as  the  prokaryotic  proto¬ 
type,  these  cells  include  a  cell  membrane,  cytoplasm,  ribosomes,  and  a  nuclear  region. 
Cytoplasm  of  prokaryotic  cells  is  the  semifluid  substance  inside  the  cell  membrane. 
Ribosomes  consist  of  ribonucleic  acid  and  protein.  Finally,  the  nuclear  region  has 
DNA  arranged  in  one  large,  circular  chromosome  [6:103]. 

The  structure  for  eukaryotic  cells  is  more  complex  than  prokaryotic  cells.  These 
type  of  cells  include  a  membrane-enclosed  cell  nucleus,  with  a  nuclear  envelope,  nu¬ 
cleoplasm,  nucleoli,  and  chromosomes  [6:103].  The  following  list  outlines  structural 
components  of  prokaryotic  cells  [6:91-92], 

•  Cytoplasm  contains  elements  of  a  cytoskeleton,  a  fibrous  network  that  give 

support  and  shape  these  cells. 
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•  Cell  nucleus  is  a  distinct  organelle  with  a  nuclear  envelope,  nucleoplasm,  nu¬ 
cleoli,  and  chromosomes. 

•  Nuclear  envelope  consists  of  a  double  membrane,  each  layer  of  which  is  struc¬ 
turally  like  the  plasma  membrane. 

•  Nucleoplasm  is  the  semifluid  portion  of  the  nucleus. 

•  Nucleoli  contain  a  significant  amount  of  RNA  and  serve  as  sites  for  the  assembly 
of  ribosomes. 

•  Chromosomes — typically  paired — contain  DNA  and  proteins  called  histones. 

•  Histones  contribute  directly  to  the  structure  of  chromosomes. 

More  information  about  intracelluar  structure  can  be  found  in  [6:103  -  104], 

2.2.3  Eukaryotic  Cell  Division  Cycle.  “To  reproduce  itself,  a  cell  must 
duplicate  all  its  components  and  separate  them,  more  or  less  evenly,  to  two  daughter 
cells,  so  that  each  daughter  has  the  information  and  machinery  necessary  to  repeat 
the  process.  In  general,  eukaryotic  cells  replicate  and  partition  their  genetic  material 
in  two  distinct,  coordinated  processes  [8:369].”  These  two  phases  are  described  below. 

S  phase  The  DNA  molecule  in  each  chromosome  is  precisely  replicated  to  form 
two  identical  sister  chromatids  that  are  held  together  by  cohesins  (tethering 
proteins). 

M  phase  The  cell  builds  a  mitotic  spindle,  condenses  its  replicated  chromosomes, 
aligns  them  on  the  midplane  of  the  spindle,  and  then,  at  anaphase,  removes 
the  cohesins  and  separates  sister  chromatids  to  opposite  poles  of  the  spindle 
[8:369], 

These  two  phases  are  separated  temporally  by  gaps  (G1  and  G2  phases)  [8:369]. 
Figure  2.1  depicts  the  eukaryotic  cell  cycle  in  general. 

Several  mechanisms  are  in  place  to  coordinate  the  proper  execution  of  events 
during  the  cell  division  cycle.  For  example,  licensing  factors  (Mcm2-7  and  Cdc6) 
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Figure  2.1.  This  diagram  of  the  typical  eukaryotic  division  cell  cycle  is  based  on  a 
figure  presented  in  [20:2],  Metaphase  is  when  chromosomes  are  aligned 
between  two  mitotic  spindles.  Anaphase  is  when  the  glue  that  holds 
the  sister  chromatids  together  is  dissolved,  allowing  each  chromatid 
to  be  pulled  by  the  microtubules  to  one  of  the  poles  of  the  spindle. 
Tclcphase  is  when  two  nuclei  are  created  just  before  the  cell  divides 
(called  cytokinesis). 
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make  sure  that  each  DNA  molecule  is  replicated  only  once  per  cycle.  Also,  cycle  logic 
does  not  allow  the  cell  to  commence  anaphase  until  DNA  replication  is  complete  and 
each  pair  of  sister  chromatids  is  properly  aligned  [8:369].  Specific  information  on  cell 
cycle  control  mechanisms  can  be  found  in  [8:369-374], 

2.2.4  Bacteriophage.  Bacteriophages  or  phages  “infect  bacterial  cells 
(hosts)  and  reproduce  within  them  [6:200]”.  Once  the  nucleic  acid  of  the  phage 
enters  the  bacterial  cell,  “further  events  follow  one  of  two  pathways,  depending  on 
whether  the  page  is  virulent  or  temperate  [6:200].”  An  infectious  virulent  phage  is 
capable  of  causing  the  death  of  a  bacteria  cell.  ’’When  the  cell  becomes  filled  with 
a  hundred  or  more  phages,  phage  enzymes  rupture  the  cell,  releasing  newly  formed 
phages,  which  can  then  infect  other  cells  [6:200].”  This  rupture  of  the  infected  cell  is 
called  a  lytic  cycle.  On  the  other  hand,  “a  temperate  phage  ordinarily  does  not  cause 
a  disruptive  infection.  Instead  the  phage  DNA  is  incorporated  into  a  bacterium’s 
DNA  and  is  replicated  with  it  [6:200].  ”  More  information  on  phages  can  be  found 
in  [6:200-203], 

2.2.5  Vibrio  fischeri.  “Vibrio  fishcheri  is  a  marine  bacterium  found  both 
as  a  free-living  organism  and  as  a  symbiont  of  some  marine  fish  and  squid  [1:8].” 
V.  fishcheri,  a  single-celled  organism,  has  a  luminescence  ability  that  is  directly 
effected  by  the  number  of  other  V.  fishcheri  cells  residing  at  a  close  proximity  (i.e., 
local  population).  As  a  free-living  organism,  population  size  is  generally  small,  and 
luminescence  appears  to  be  absent.  However,  as  a  symbiont,  population  size  is 
generally  large  (dense)  and,  luminescence  is  usually  detectable. 

Luminescence  in  V.  fischeri  is  activated  by  its  quorum  sensing  system.  When 
the  local  popluation  reaches  a  quorum — a  minimum  population  size — the  quorum 
sensing  system  activates  a  set  of  genes  that  enable  the  cell  to  become  lumines¬ 
cent  [1:8].  For  information  concerning  genes,  see  section  2.3.2.  More  information 
about  V.  fishcheri  can  be  found  in  [1:7-8]. 
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2.3  Biochemistry 

Biochemistry  is  the  study  of  chemical  reactions  that  occur  in  living  systems. 
Biochemistry  focuses  specifically  on  the  molecules  of  matter  (biomolecules),  which 
living  organisms  are  composed  of.  Four  main  classes  of  important  biomolecules  are 
proteins,  carbohydrates,  lipids,  and  nucleic  acids  [27:1].  Key  objects  in  the  modeling 
process  are  encompassed  in  this  held  of  study. 

Dynamic  behavior  of  the  intracellular  environment  is  characterized  by  flux 
of  various  types  of  molecules  (species)  that  are  changed  by  biochemical  reactions. 
Understanding  the  principles  of  biochemistry  is  essential  in  being  able  to  incorporate 
this  dynamic  behavior  into  candidate  models.  Some  key  principles  of  biochemistry 
are  the  functionality  of  proteins  and  nucleic  acids,  how  biochemical  reactions  take 
place,  and  how  chemical  kinetics  govern  reactions.  I  discuss  these  topics  in  the 
following  sections. 

2.3.1  Proteins.  In  the  book,  Evolutionary  Computation  in  Bioinformatics, 
the  author  describes  the  significance  of  protein  molecules,  “Proteins  make  up  most 
of  an  organism’s  biomass,  and  play  a  key  role  in  its  metabolic  and  other  cellular  and 
bodily  processes  [9:12].”  Specifically,  a  protein  is  a  linear  chain  of  amino  acids.  The 
amino-acid  sequence  of  the  protein,  along  with  its  three-dimensional  shape,  is  the 
functional  portion  of  information  flow  in  a  cell  [14:545].  An  enzyme  is  a  special  type 
of  protein  that  is  crucial  in  enabling  chemical  reactions  to  take  place. 

An  enzyme  is  a  catalyst  for  chemical  reactions  in  living  organisms  [27:9].  Dur¬ 
ing  chemical  reactions,  a  specific  enzyme  associates  with  a  corresponding  substrate — 
forms  a  complex — changes  it  into  a  product,  and  then  dissociates  with  the  new 
product.  However,  the  enzyme  remains  unchanged  and  is  ‘free’  to  take  part  in  an¬ 
other  chemical  reaction  with  a  specific  substrate.  Enzyme  specificity  is  related  to 
its  shape.  Detailed  information  on  proteins  and  enzymes  can  be  found  in  [6:41-42] 
and  [27:9-18], 
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2.3.2  Nucleic  Acids.  Nucleic  acids  are  the  carriers  of  the  genetic  code. 
“They  contain  genetic  information  that  determines  all  the  heritable  characteristics 
of  a  living  organism,  be  it  a  microbe  or  a  human.  Such  information  is  passed  from 
generation  to  generation  and  directs  the  protein  synthesis  in  each  organism  [6:43].” 
RNA  and  DNA  are  the  two  nucleic  acids  found  in  living  organisms  [6:44],  The  data 
structure  for  genetic  information  storage  within  each  cell  is  stored  in  the  following 
hierarchical  configuration — genome,  chromosomes,  and  genes. 

The  genome  contains  the  complete  set  of  instructions  for  an  organism.  “It 
contains  the  master  blueprint  for  all  cellular  structures  and  activities  for  the  lifetime 
of  the  cell  or  organism  [13:5].”  The  human  genome  is  subdivided  into  structures 
called  chromosomes. 

Chromosomes  consist  of  tightly  coiled  threads  of  DNA  and  associated  protein 
molecules.  Apart  from  reproductive  cells  and  mature  red  blood  cells,  every  cell  in 
the  human  body  contains  23  pairs  of  chromosomes.  Each  chromosome  is  a  packet 
of  compressed  and  entwined  DNA.  The  human  genome  contains  three  billion  base 
pairs  [13:5].  Subunits  of  a  chromosomes  are  genes. 

Genes  are  the  basic  physical  and  functional  units  of  heredity.  They  are  charac¬ 
terized  as  a  “specific  sequence  of  nucleotide  bases,  whose  sequences  carry  the  infor¬ 
mation  required  for  constructing  proteins,  which  provide  the  structural  components 
of  cells  and  tissues  as  well  as  enzymes  for  essential  biochemical  reactions  [13:7].”  The 
human  genome  contains  about  30,000  to  40,000  genes  [1:2].  Specific  information  on 
nucleic  acids  can  be  found  in  [13:5-9]  and  [6:42-45]. 

2.3.3  Chemical  Reactions.  Chemical  reactions  are  continuously  taking 
place  within  a  cell.  As  a  result,  a  cell  is  like  a  factory,  where  products  are  continuously 
being  produced.  Of  course,  this  production  is  highly  coordinated.  As  indicated 
in  the  following  quote,  enzyme  regulation  is  a  main  factor  in  regulating  chemical 
reactions  [6:109]. 


2-9 


Like  nearly  all  other  chemical  processes  in  living  organisms... each  consist 
of  a  series  of  chemical  reactions  in  which  the  product  of  one  reaction 
serves  as  the  substrate  (reacting  material)  for  the  next.  Each  reaction  in 
a  pathway  is  controlled  by  a  particular  enzyme. 

Hasty  et  al.  discusses  the  regulation  of  intracellular  chemical  reactions  from  a 
standpoint  of  gene  regulatory  networks.  They  explicitly  define  negative  and  positive 
feedback  as  regulatory  modes  within  a  cellular  network.  Negative  feedback  is  when 
a  cellular  network  inhibits  in  own  level  of  activity.  Contrarily,  positive  feedback  is 
when  such  a  network  increases  its  own  level  of  activity  [10:269-270].  The  models 
of  Vibrio  hscheri  (a  unicellular  bacteria)  and  saccharoyces  cerevisiae  (yeast  cell)  are 
covered  in  this  thesis.  They  explicitly  use  positive  and  negative  feedback  to  regulate 
the  production  of  end  products. 

2.3.4  Chemical  Kinetics.  Intracellular  chemical  reactions  occur  at  various 
rates  involving  a  wide  diversity  of  molecular  interactions.  Consequently,  RNA  and 
proteins  are  in  flux.  In  order  to  build  meaningful  intracellular  network  models,  this 
flux  must  be  captured  in  candidate  models.  Thus,  a  fundamental  understanding 
of  chemical  kinetics  is  necessary.  In  this  section,  some  chemical  kinetic  terms — 
chemical  reaction  equation  notation,  reaction  rates,  and  reaction  rate  laws-are  briefly 
discussed. 

A  chemical  reaction  can  be  assigned  a  specific  reaction  rate  (velocity)  at  which 
that  chemical  reaction  takes  place.  Reaction  rate  is  defined  “as  the  change  in  con¬ 
centration  per  unit  of  time  of  any  reaction  product  for  which  the  stoichiometric  co¬ 
efficient  is  1  (assuming  the  volume  of  the  reacting  system  remains  constant)  [25:4].” 
Stoichiometry  is  “the  determination  of  the  proportions  in  which  chemical  elements 
combine  or  are  produced  and  the  weight  relations  in  any  chemical  reaction  [19:1320].” 

Various  units  of  measure  are  used  to  quantify  concentration  in  chemical  re¬ 
actions.  Some  common  units  of  measure  include  moles  per  liter  moles  per 
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cubic  centimeter  (7^),  and  molecules  per  cubic  centimeter  ( ",ol^les )  [25:4],  Plus, 
the  second  is  generally  the  time  unit. 


For  example,  I  generalize  an  example  taken  from  [25:5].  The  following  one-way 
reaction  in  some  unit  of  measure  is  assumed  at  a  constant  volume  and  obtain  a  rate 
of  change  of  B2  concentration  of  ■ 


2 AB  ->  Ao  +  Bo 


(2.1) 


Then  that  is  the  rate  of  the  reaction.  From  stoichiometry,  the  other  reaction  rates 
can  be  deduced. 

d[A2]  1.5  mole 


and 


dt 


d[AB2 


10 5lsec 

—3.0  mole 


(2.2) 


(2.3) 


dt  10  Hsec 

In  general,  “for  any  reaction  product,  the  rate  of  change  of  concentration  is  equal 
to  the  stoichiometric  coefficient  times  the  rate  of  reaction;  while  for  the  reactants, 
the  rate  of  change  of  concentration  is  equal  to  the  negative  of  the  stoichiometric 
coefficient  times  the  rate  of  reaction  [25:5].” 

Furthermore,  the  first  equation  can  be  rearranged  so  that  the  substances  are 
all  on  the  left  hand  side,  with  negative  stoichiometric  coefficients  for  the  reactions. 


A2  +  B2  -  2 AB  =  0. 


(2.4) 


In  the  generalized  form,  this  equation  can  be  written  as 


E  Mi  =  0,  (2.5) 

where  Ai  is  the  symbol  of  a  product  or  reactant,  and  vt  is  its  stoichiometric  coeffi¬ 
cient.  Plus,  Vi  is  positive  for  products  and  negative  for  reactants.  Now,  with  this 
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redefinition  of  the  stoichiometric  coefficients  and  R  representing  the  rate  of  reaction, 
then 

IT  =  *'<*•  <26) 

Equations  (1.2)  and  (1.3)  can  be  derived  using  equation  (1.4).  In  this  case,  both 
stoichiometric  coefficients  equal  one  and  R  =  • 

When  building  intracellular  network  models,  a  specific  reaction  rate  law  is  often 
assumed  for  each  reaction.  At  a  given  temperature,  and  perhaps  within  a  limited 
range  of  concentrations,  one  can  write  a  rate  law  for  a  reaction  of  the  form, 


JtnwW  <2-7> 

i 

where  the  Aj  are  the  reactants,  the  X j  are  substances  that  are  not  reactants  but  do 
influence  the  rate,  the  a’s  and  /3’s  are  coefficients,  that  are  not  necessarily  related 
to  the  stoichiometric  coefficients  u,  and  k  is  a  rate  constant  [25:9].  However,  in  this 
thesis  I  will  only  consider  elementary  reactions.  Therefore,  “for  many  elementary 
reactions,  the  rate  law  coefficients  are  equal  to  the  stoichiometric  coefficients  [25:9].” 
Specifically,  only  mass  action  and  Michaelis-Menten  rate  laws  are  implemented  in 
this  thesis. 

Mass  action  rate  laws  are  expressed  in  the  form  of  equation  2.7.  However,  for 
mass  action  rate  laws,  a’s  and  /3’s  are  stoichiometric  coefficients  for  their  respective 
substances.  Unlike  this  rate  law,  the  Michaelis-Menten  kinetic  rate  law  explicitly 
accounts  for  reactions  catalyzed  by  enzymes. 

“Most  biological  reactions  are  modeled  by  the  Michaelis-Menten  kinetic 
scheme  [9:258].”  Michaelis-Menten  rate  law  is  characterized  by 

h[S][E]0 
[5]  +  Km 
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where  S  is  the  reactant  or  substrate  molecule,  [E] o  is  the  amount  of  enzymatic 
concentration,  Km  is  the  Michaelis  constant,  and  k 2  is  the  rate  constant  for  the 
combination  of  an  enzyme  and  an  (unreacted)  substrate  molecule  ( ES )  becoming 
dissociated  [25:131]. 

Closely  following  [25:130-131],  the  Michaelis-Menten  equation  is  derived  with 
the  assumption  that  the  concentration  of  \E\  is  significantly  smaller  than  [S']  (i.e., 
[E]  <C  [S']).  Since  [E]  <C  [S'],  a  steady  state  concentration  of  [.ES]  will  quickly  build 
up.  Therefore,  rate  for  this  reaction  is  V  =  A^fES].  This  equation  depicts  the 
appearance  of  the  product  P  as  E  catalyzes  S  (discussed  below). 

The  object  in  deriving  the  Michaelis-Menten  equation  is  to  represent  V  in 
terms  of  experimentally  measurable  quantities  of  [S']  and  [E]0.  This  is  accomplished 
in  the  following  derivation. 

The  author  starts  with  a  simple  unimolccular  reaction, 

S  ->  P  (2.9) 

where  S  is  the  substrate  molecule,  and  P  is  the  product  of  the  reaction.  This  reaction 
can  then  be  divided  into  two  different  stages,  represented  by  the  following  pair  of 
reaction  equations. 

fci 

E  +  Sr.ES  (2.10) 

where  k\  is  the  forward  rate  constant  and  k-\  is  the  reverse  reaction  rate  constant. 

ES^E  +  P  (2.11) 

where  is  the  forward  reaction  rate  constant. 

In  equation  2.10  (moving  left  to  right),  S'  reversibly  associates  with  E.  In 
equation  2.11 — the  enzyme-substrate  complex  ES — breaks  down  into  the  original 
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enzyme  and  a  product.  E  can  be  used  over  and  over  again  to  catalyze  additional  S 
molecules.  This  equation  represents  appearance  of  P  after  E  catalyzes  S. 

If  equation  2.9  is  a  two-way  reaction  instead  of  a  one-way  reaction,  the  following 
can  be  deduced  from  equations  2.9,  2.10,  and  2.11. 

S^P  (2.12) 

Again,  this  reaction  is  divided  into  two  different  stages,  represented  by  the  following 
pair  of  reaction  equations. 

ki 

E  +  SlP.ES  (2.13) 

where  k\  is  the  forward  reaction  rate  constant  and  k_\  is  the  reverse  reaction  rate 
constant. 

fc2 

ESkP2E  +  P  (2.14) 

where  k2  is  the  forward  reaction  rate  constant  and  /c_2  is  the  reverse  reaction  rate 
constant. 

This  new  set  of  equations  is  the  same  as  the  above  set  with  the  exception  of 
2.14.  In  this  case,  equation  2.14  is  reversible.  As  a  result,  enzyme  and  a  product — E 
and  P — can  bind  back  together  to  again  form  the  complex  ES. 

Using  equations  2.10  and  2.11  and  applying  standard  techniques  of  biochemical 
kinetics  ([25:8-10]),  the  following  reaction  velocities  can  be  expressed, 

v1  =  -ki  [£][S]  +  k-i[ES]  (2.15) 

v2  =  - k2[ES }  (2.16) 
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where  each  term  in  equation  2.15  is  obtained  using  equation  2.7.  Therefore,  the  rate 
of  change  for  [ES]  can  be  expressed  as  follows. 


d[ES] 

—7—  =  -V!  +  V2 

at 

=  ki [E] [S']  -  h.ilES]  -  k2[ES }  (2.17) 

There  is  also  a  conservation  relationship  between  the  total  concentration  of 
enzyme  [77]  0  and  the  sum  of  the  concentration  of  enzyme  associated  with  the 
substrate  [77  S']  and  free  enzyme  [E],  This  relationship  is  expressed  as 

[E] o  =  [ES]  +  [E]  =>  [E]  =  [E] o  -  [ES]  (2.18) 

Since  [E]  [S'],  this  implies  a  steady  state  concentration  of  [Sis']  will  quickly  occur 

[25:130].  This  means  that  the  rate  of  the  overall  reaction  is 

V  =  k2[ES[.  (2.19) 

Therefore,  [ES]  remains  constant,  which  implies  =  0.  Therefore, 


0  =  ki[E][S]  -  k-i[ES]  -  k2[ES ] 


(2.20) 


[ES](k_1  +  k2)  =  k1[E][S] 


[ES] 


hjS] 

(k-i  +  k2 ) 


[E] 
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[ES](k_1  +  k2)  =  k1[S}([E\0-[ES]) 


[ES](k.1  +  k2  +  k1[S})=k1[S][E}o 


[Po1  MgHgk 
[  J  (k-1  +  k2  +  k1[s]) 

Substituting  this  into  equation  2.19,  in  which  the  velocity  of  the  build  up  of  [ ES ]  is 
expressed,  reduces  to  the  following  equation. 

,,  Mgligk.  Mdgk  ,,,,, 

2  (k-i  +  k2  +  ki[S])  [S]  +  (k-i  +  k2) / ki 


The  author  then  expresses  the  second  term  in  the  denominator  as  Krn  = 
(k- 1  +  k2)/k\ — the  Michaelis  constant.  Substituting  Km  into  equation  2.21  com¬ 
pletes  the  derivation  of  the  Michaelis-Menten  equation.  Additionally,  the  author 
discusses  limiting  behavior  at  low  and  high  substrate  concentrations  exhibited  by 
the  Michaelis-Menten  equation  [25:131].  At  low  values  of  [S'], 


V 


k2[S}[E]  o 


Kn 


(2.22) 


which  gives  a  first-order  dependence  on  [S'].  In  this  case,  there  are  a  lot  of  ‘free’ 
enzyme  molecules.  As  [S']  is  increased,  the  rate  of  reaction  is  increased,  up  until  all 
enzyme  molecules  are  associated  with  substrate  molecules  [25:131]. 
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On  the  other  hand,  at  high  levels  of  [S'],  this  term  rather  than  Km  dominates 
the  denominator  to  give 

V  =  h[E]0  (2.23) 

In  this  case,  the  rate  of  reaction  no  longer  depends  on  [S].  When  this  state  occurs, 
all  enzymes  are  associated  with  substrate  molecules  (no  free  enzymes).  Hence,  the 
reaction  has  reached  its  maximum  velocity  ( Vmax )•  Additional  substrate  molecules 
will  not  be  catalyzed  until  enzymatic  molecules  disassociate  with  substrate  molecules 
currently  undergoing  chemical  change  [25:131]. 

This  Michaelis-Menten  equation  is  implemented  in  chapters  three  and  four. 
Using  the  rate  equation  approach  (discussed  later),  it  is  expressed  in  several  reaction 
equations  of  the  budding  yeast  cell  cycle  model. 

2.4  Biological  Pathways 

A  metabolic  pathway  consists  of  a  series  of  chemical  reactions  in  which  the 
product  of  one  reaction  serves  as  the  substrate  (reacting  material)  for  the  next 
chemical  reaction.  Each  reaction  in  a  metabolic  pathway  is  controlled  by  a  par¬ 
ticular  enzyme.  For  example,  A^B^C^D-^E  represents  a  simple  metabolic 
pathway.  In  this  pathway,  A  is  the  initial  substrate,  E  is  the  final  product,  and  B, 
C ,  and  D  are  intermediates  [6:109]. 

Metabolic  pathways  are  either  catabolic  or  anabolic  (biosynthetic).  Catabolic 
pathways  capture  energy  in  a  form  cells  can  use.  Anabolic  pathways  make  the  com¬ 
plex  molecules  that  form  the  structure  of  cells,  enzymes,  and  other  molecules  that 
control  cells.  These  pathways  use  building  blocks  such  as  sugars,  glycerol,  fatty 
acids,  amino  acids,  nucleotides,  and  other  molecules  to  make  carbohydrates,  lipids, 
proteins,  nucleic  acids,  or  a  combination  of  these  molecules.  Adenosine  triphosphate 
(ATP)  molecules  are  the  links  that  couple  catabolic  and  anabolic  pathways.  Energy 
released  in  catabolic  reactions  is  captured  and  stored  in  the  form  of  ATP  molecules, 
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which  arc  later  broken  down  to  provide  the  energy  needed  to  build  up  new  molecules 
in  biosynthetic  pathways  [6:109].  Glycolysis  is  an  example  of  a  metabolic  path¬ 
way  [6:114],  The  glycolytic  pathway  model  in  section  A  incorporates  ATP  molecules. 

2. 5  Phosphorylation 

“Phosphorylation  is  the  addition  of  a  phosphate  group  to  a  molecule,  often 
from  ATP.  This  addition  generally  increases  the  molecule’s  energy.  Thus,  phosphate 
groups  commonly  serve  as  energy  carriers  in  biochemical  reactions  [6:116].'"  In  the 
glycolytic  pathway  model  in  section  A,  phosphate  groups  from  two  molecules  of  ATP 
are  added  to  glucose  at  the  beginning  part  of  the  pathway.  This  expenditure  of  two 
ATPs  raises  the  energy  level  of  glucose,  enabling  it  to  then  participate  in  subsequent 
reactions  [6:116]. 

Adenosine  diphosphate  (ADP)  has  a  direct  correlation  with  ATP  in  regards  to 
glycolysis.  With  ADP  and  inorganic  phosphate  (Pi)  available  in  the  cytoplasm,  the 
energy  released  from  substrate  molecules  is  used  to  form  high-energy  bonds  between 
ADP  and  Pf 

ADP  +  Pi  +  energy  ATP  (2.24) 

In  this  way  energy  is  captured  in  ATP  at  the  substrate  level  [6:116].  Phosphoryla¬ 
tion  is  a  key  activity  in  the  glycolytic  pathway  and  in  the  budding  yeast  cell  cycle 
(Saccharomyces  cerevisiae). 

2.6  Protein  Synthesis 

With  an  understanding  of  the  fundamental  principles  of  Biochemistry,  it  is  fea¬ 
sible  to  analyze  one  of  the  most  important  activities  of  intracelluar  function — protein 
synthesis.  As  previously  mentioned,  proteins  are  crucial  components  in  major  cell 
functions.  The  degree  of  importance  is  quantified  by  [6:169]:  “All  cells  must  con¬ 
stantly  synthesize  proteins  to  carry  out  their  life  processes:  reproduction,  growth, 
repair,  and  regulation.”  Protein  synthesis  involves  several  orchestrated  steps:  the 
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gene  (or  groups  of  genes)  is  tagged  for  transcription  (initialized),  the  gene  is  tran¬ 
scribed  to  a  messenger  RNA  (mRNA),  and  mRNA  interacts  with  a  ribosome  in 
building  a  specified  protein  from  a  sequence  of  amino  acids  previously  prescribed  by 
the  gene  [9:11].  In  this  manner,  a  gene  or  groups  of  genes  indirectly  build  specified 
proteins. 

Protein  synthesis  is  initiated  when  “a  pre- initiation  complex  is  assembled  around 
the  promoter  region  just  upstream  of  the  gene  that  is  to  be  expressed  [9:11].”  This 
complex  attracts  RNA  polymerase,  a  protein  that  causes  the  strands  of  DNA  to 
separate  near  the  start  of  the  gene.  Moving  along  this  gene  or  exon  (the  active  part 
of  DNA  transcribing  mRNA)  RNA  polymerase  affects  production  of  an  RNA  tran¬ 
script  [9:11].  Next,  this  mRNA  transcript  interacts  with  a  ribosome — the  location 
where  protein  is  synthesized  via  translation  of  mRNA.  In  prokaryotic  cells,  tran¬ 
scription  and  translation  takes  place  both  in  the  cytoplasm,  whereas  in  eukaryotes, 
transcription  takes  place  within  the  cell  nucleus,  while  translation  takes  place  outside 
the  nucleus  [6:171].  Based  on  [6:176],  the  following  steps  summarize  translation. 

1.  An  mRNA  transcript  becomes  properly  oriented  on  a  ribosome. 

2.  A  ribosome  reads  each  codon  of  the  mRNA,  and  the  appropriate  tRNA  com¬ 
bines  with  it  and  delivers  a  particular  amino  acid  to  the  ribosome  (a  codon  is  a 
sequence  of  three  bases  in  mRNA  that  specifies  a  particular  amino  acid  [6:G13]). 

3.  Step  2  is  repeated  until  a  stop  codon  is  read  that  terminates  protein  synthesis. 

Using  a  diagram  in  [9:13]  as  a  reference,  Figure  2.2  illustrates  how  a  protein  is 
generated  from  genetic  information. 

2. 7  Mathematical  Modeling 

“Mathematical  models  are  useful  for  providing  a  framework  for  integrating 
data  and  gaining  insights  into  the  static  and  dynamic  behavior  of  complex  biological 
systems  [26:247]....”  This  class  of  models  approximates  the  dynamics  of  biological 
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Promoter 


Gene 


F  olds  into  a  protein 


Figure  2.2.  Steps  (a)-(f)  illustrate  how  a  protein  is  made  from  genetic  information 
(transcription/translation),  (a)  Section  of  DNA  containing  a  gene, 
preceded  by  a  promoter  region,  (b)  An  RNA  polymerase  complex  is 
attracted  to  the  DNA  around  the  promoter  and  begins  to  interact  with 
the  DNA  at  the  start  of  the  gene-beginning  of  transcription,  (c)  The 
RNA  ploymerase  moves  along  from  left  to  right,  gradually  building  an 
RNA  transcript  for  the  gene  sequence,  (d)  The  transcript  is  complete, 
and  it  now  separates  from  both  the  RNA  polymerase  and  the  DNA.  (e) 
The  ribosome  uses  it  to  automatically  manufacture  a  protein  molecule- 
translation  in  progress,  (f)  Translation  is  complete.  This  diagram  is 
based  on  an  illustration  published  in  [9:13]. 
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systems  using  a  system  of  ODEs  (typically  nonlinear)  as  the  computational  model. 
In  addition,  the  system  of  ODEs  can  be  augmented  with  logic  terms  that  capture 
events  (e.g.,  cell  division).  Hasty  et  al.  defines  the  following  states  and  regulative 
mechanisms  of  biological  networks  that  have  been  thoroughly  studied  [10:269-270]. 
It  is  necessary  to  have  an  understanding  of  these  terms,  because  these  traits  and 
regulative  mechanisms  are  normally  incorporated  into  candidate  models. 

Equilibrium  state:  For  a  gene  product  (such  as  a  protein),  the  rate  of  production 
is  balanced  by  its  rate  of  degradation  [10:269]. 

Fixed  point:  A  point  at  which  the  rates  of  change  of  all  variables  in  a  system  are 
exactly  zero.  A  system  precisely  at  its  fixed  point  (or  steady  state)  will  remain 
there  permanently  [10:269]. 

Multistability:  The  property  of  having  more  than  one  stable  fixed  point  [10:269]. 

Negative  feedback:  A  component  (species)  of  a  system  is  subject  to  negative  feed¬ 
back  when  it  inhibits  its  own  level  of  activity  [10:269]. 

Positive  feedback:  A  component  (species)  of  a  system  is  subject  to  positive  feed¬ 
back  when  it  increases  its  own  level  of  activity  [10:270]. 

Certain  configurations  of  negative  and  positive  feedback  can  lead  to  system 
multistability.  This  is  an  important  concept  for  the  modeler  to  grasp,  because  “un¬ 
derstanding  how  multistability  arises  is  thus  relevant  to  understanding  the  opera¬ 
tion  of  natural  biological  switches,  as  well  as  the  design  of  synthetic  switching  net¬ 
works  [10:269].”  Specific  details  on  multistability  and  its  implications  can  be  found 
in  [10:269]. 

2.8  Modeling  Approaches 

Once  the  modeler  has  obtained  a  well-rounded  understanding  on  the  basics  of 
cellular  structure,  biochemistry,  and  mathematical  modeling  concepts,  the  modeler  is 
nearly  ready  to  consider  a  modeling  approach  for  implementation.  Prior  to  selecting 
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a  modeling  approach,  the  modeler  should  be  familiar  with  modeling  alternatives.  I 
discuss  three  alternative  modeling  approaches:  Boolean,  chemical  kinetics  or  rate 
equation,  and  stochastic  [10:270].  These  approaches  range  from  low  to  high  level  of 
detail  in  characterizing  a  biological  network. 

Starting  from  a  low  level  of  network  detail,  the  Boolean  approach  enables  con¬ 
trol  of  network  dynamics  by  assigning  a  switch  to  each  gene.  In  this  way,  each  gene 
can  assume  only  one  of  two  states — ON  or  OFF.  Also,  genes  within  the  network  have 
the  ability  to  interact  with  each  other  through  user-defined  interaction  rules  [10:270]. 

For  example,  suppose  three  genes — A,  B,  and  C — collectively  regulate  a  net¬ 
work.  A  possible  rule  for  network  dynamics  could  be:  if  genes  A  and  B  are  OFF 
at  the  preceding  time  step,  then  gene  C  is  turned  ON  during  the  current  time  step. 
This  represents  synchronous  control  of  genetic  states.  Asynchronous  methods  for 
activating  and  deactivating  gene  switches  can  also  be  incorporated. 


The  advantage  of  the  Boolean  approach  is  that  such  models  are  not  difficult 
to  implement.  Also,  the  computational  cost  to  run  simulations  is  relatively  low. 
The  disadvantage  of  this  approach  is  that  “the  abstraction  of  genes  to  ON/OFF 
switches  makes  it  difficult  or  impossible  to  include  many  of  the  details  in  cellular 
biology  [10:270]”. 


The  chemical  kinetics  or  rate  equation  approach  approximates  a  cellular  net¬ 
work  at  a  midrange  level  of  detail.  As  stated  in  section  2.7,  a  system  of  ODEs  (typ¬ 
ically  nonlinear)  approximates  network  dynamics.  In  this  case,  real-valued  species 
concentration  variables  are  updated  by  approximating  a  solution  vector  in  respect  to 
the  current  timestep.  In  a  simple  network,  the  rate  of  change  for  each  species  concen¬ 
tration  is  equated  to  concentration  levels  of  other  interacting  species  concentrations. 
Each  interaction  has  a  positive  or  negative  effect  on  the  rate  of  concentration  change 
for  a  given  species.  Equation  2.25  mathematically  defines  this  interaction  among 
species. 


dXj 

dt 


v;{x)  -  V&X) 


(2.25) 
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where  X*  is  the  vector  of  species  concentrations,  V*  and  VJ  are  the  kinetic  rate  raws 
previously  mentioned  in  the  section  2.3.4  that  produce  and  consume  Xt  respectively. 

The  main  advantage  of  the  rate  equation  approach  is  that  dynamics  are  totally 
characterized  by  a  system  of  ODEs.  Accordingly,  the  modeler  can  readily  analyze 
and  predict  model  behavior  by  applying  numerical  analysis  methods  to  the  system 
of  ODEs.  For  example,  stability  of  steady  state  points  can  be  determined  by  cal¬ 
culating  the  eigenvalues  of  the  Jacobian  matrix  of  the  nonlinear  system.  The  main 
disadvantage  for  this  method  is  that  it  is  totally  deterministic.  Consequently,  this 
approach  ignores  the  inherent  randomness  or  ‘noise’  that  acts  upon  biological  sys¬ 
tems.  However,  ‘noise’  can  be  added  to  this  deterministic  model  by  augmenting  it 
with  stochastic  terms. 

Unlike  the  deterministic  rate  equation  approach,  the  stochastic  kinetics  ap¬ 
proach  that  incorporates  the  highest  level  of  detail  explicitly  takes  into  account  the 
randomness  or  fluctuations  in  rates  of  gene  expression.  Results  from  Arkin  et  al. 
suggest  that  such  random  behavior  “can  produce  a  highly  erratic  time  pattern  of 
protein  production  in  each  individual  cell  and  a  wide  diversity  of  protein  concen¬ 
trations  across  a  cell  population  at  any  instant  in  time  [3:1633].”  Arkin  et  al. 
state  the  importance  of  considering  stochastic  behavior:  “When  two  independently 
produced  regulatory  proteins  acting  at  low  cellular  concentrations  competitively  con¬ 
trol  a  switch  point  in  a  pathway,  stochastic  variations  in  their  concentrations  can 
produce  probabilistic  pathway  selection....  [3:1633]”  In  short,  the  stochastic  kinetics 
approach  attempts  to  closely  model  important  cellular  events  that  are  inherently 
random  within  a  genetic  regulatory  network.  This  method  enables  the  modeler  to 
track  events  at  the  micro  level.  Detailed  information  on  the  stochastics  approach 
can  be  found  in  [3:1633-1648]. 

Likewise,  this  approach  has  a  significant  advantage  and  major  disadvantage. 
“[It]  is  impressively  complete  and  yields  a  detailed  picture  of  the  behavior  of  the 
system  modeled.”  The  major  disadvantage  is  that  model  results  “comes  at  a  high 
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computational  cost  and  sacrifices  any  immediate  prospect  of  analytical  treatment 
[10:270].” 

2.9  Gepasi 

“Gepasi  (version  3.21)  is  a  software  package  for  biochemical  systems.  It  simu¬ 
lates  the  kinetics  of  systems  of  biochemical  reactions  and  provides  a  number  of  tools 
to  fit  models  to  data,  optimize  any  function  of  the  model,  peform  metabolic  control 
analysis,  and  linear  stability  analysis  [17:1]”.  Gepasi  uses  the  chemical- reaction- 
centered  approach  in  building  the  computational  model.  In  this  mode,  the  user 
enters  a  system  of  reaction  equations  using  standard  chemical  notation.  From  this 
user  input,  the  program  automatically  builds  and  solves  the  system  of  ODEs.  Unfor¬ 
tunately,  the  user  is  unable  to  view  the  system  of  ODEs  generated  by  the  program. 

After  building  the  model,  the  user  can  perform  two  kinds  of  simulations — time 
course  and  steady  state.  Time  course  simulations  depict  how  kinetics  of  the  reaction 
network  evolve  over  time  from  a  starting  point.  In  setting  up  this  type  of  simulation, 
the  user  defines  the  end  time,  number  of  sampling  points,  and  output  configuration 
for  the  resultant  data  hie.  Output  hies  are  configured  in  such  a  way  that  third  party 
applications  can  readily  import  Gepasi  data  hies  [17:1]. 

Steady  state  simulations  attempt  to  find  a  steady  state  near  a  starting  point. 
Gepasi  cannot  find  some  types  of  steady  states  (e.g.,  saddle  points).  Steady  state  re¬ 
sults  are  sent  to  a  data  hie.  These  hies  contain  calculated  steady  state  concentration 
for  each  metabolite  via  a  row  vector  [17]. 

Additional  features  listed  in  the  Gepasi  user  manual  are  ‘scanning  parameter 
space’  and  ‘fitting  model  parameters  to  experimental  data’.  In  the  scanning  pa¬ 
rameter  mode,  the  program  automatically  performs  a  sensitivity  analysis  on  model 
parameters  as  specihed  by  the  user.  In  the  fitting  mode,  the  program  minimizes  the 
sum  of  squares  of  residuals  between  predicted  values  from  the  model  and  experimen¬ 
tal  data  [17]. 
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2.10  Virtual  Cell 


Virtual  Cell  is  another  software  package  for  modeling  cellular  biological  pro¬ 
cesses.  “It  is  based  on  the  mapping  of  experimental  biochemical  and  electrophys- 
iological  data  onto  experimental  images.  The  framework  is  designed  to  enable  the 
construction  of  complex  general  models  that  encompass  the  general  class  of  prob¬ 
lems  coupling  reaction  and  diffusion  [23:228].”  Unlike  Gepasi,  a  graphical  approach 
is  incorporated  into  Virtual  Cell.  In  this  approach,  the  user  constructs  a  model 
by  placing  and  manipulating  abstract  modeling  objects  onto  an  active  workspace. 
Modeling  objects  can  be  edited,  viewed,  stored  in  a  remote  database,  and  analyzed 
using  the  WWW-based  user  interface  [23:230]. 

For  building  models  and  running  simulations,  the  user  interacts  with  the  Vir¬ 
tual  Cell  application  through  a  GUI  enabling  access  to  various  editors.  After  running 
the  simulation,  the  user  can  view  results  (in  graphical  and  tabular  format)  in  separate 
windows.  Program  output  includes  ODEs,  graphical  display  of  simulation  results, 
and  simulation  data  in  tabular  format.  The  user  has  several  choices  of  formats  in  ex¬ 
porting  the  data.  Additional  information  about  Virtual  Cell  and  access  to  its  WWW- 
based  user  interface  is  at  http://www.nrcam.uchc.edu/vcellR3/login/tutorial.jsp. 

2.11  E-CELL2 

E-CELL  is  a  software  package  that  enables  the  user  to  model  metabolic  path¬ 
ways  and  higher-order  cellular  processes,  such  as  protein  synthesis  and  signal  trans¬ 
duction  [28].  This  package  employs  the  ‘Substance- Reactor  Model’  in  depicting  the 
structure  of  a  cell  and  the  chemical  reactions  that  take  place  within  a  cell.  In  this 
fashion,  three  classes  of  objects  are  used  to  model  cell  activity:  substances,  reactors, 
and  systems.  Reactors  calculate  the  changes  in  the  amount  of  substances  over  time. 
Systems  express  the  location  of  reactors  and  substances.  Reactors  implement  the 
specified  reaction  mechanism  (e.g.,  Michaelis-Menten  rate  law)  [15]. 
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The  user  can  view  simulation  results  within  several  windows.  The  state  of  a 
simulation  run  can  be  captured  and  subsequently  saved.  The  user  is  able  to  save 
simulation  results  and  whole  cell  state  information  of  the  model  at  a  given  time 
to  a  hie.  Contents  in  the  data  hie  include  simulation  data  on  all  substances  and 
reactors,  time,  quantity,  an  average  quantity,  the  maximum  quantity,  the  minimum 
quantity,  the  average  concentration,  the  maximum  concentration,  and  the  minimum 
concentration  [15].  Additional  information  can  be  accessed  from  [28]  and  [15]. 

2.12  BioSPICE 

The  Defense  Advanced  Research  Projects  Agency’s  (DARPA)  Bio-Computation 
program  is  aimed  at  exploring  and  developing  computational  methods  and  models 
at  the  biomolecular  and  cellular  levels.  This  organization  released  a  software  bundle 
called  BioSPICE  (version  1.1.).  It  contains  an  assortment  of  biochemical  mod¬ 
eling  tools.  Software  packages  included  in  this  bundle  are — among  others — Grass, 
Simpathica,  Pathway  Builder,  Charon,  Biosketchpad,  and  JigCell. 

Charon  and  Biosketchpad  are  used  in  conjunction  with  each  other,  forming 
BioCharon.  This  application  incorporates  a  graphical  approach.  In  contrast,  Jigccll 
incorporates  a  spread  sheet  approach.  I  describe  BioCharon  and  JigCell  in  chapter 

III.  I  list  and  analyze  simulation  results  produced  by  these  software  tools  in  chapter 

IV. 
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III.  Intracellular  Modeling 


3. 1  Overview 

The  intracellular  environment  involves  intricate  dynamics.  Understanding  the 
underlying  logic  that  controls  these  dynamics  is  the  foundation  for  illuminating  over¬ 
all  cell  functionality.  Despite  the  emergence  of  powerful  techniques  in  sequencing 
genes  and  profiling  RNA  and  protein  at  the  genomic  level,  “there  is  a  considerable 
gap  between  the  availability  of  new  seqence  data  and  a  scientific  understanding  of 
that  information  [9:3].”  Mathematical  modeling  can  help  bridge  that  gap  between 
the  large  repository  of  data  and  the  scientific  understanding  of  that  information. 

As  previously  mentioned  in  section  2.1,  formulating  a  computational  model 
involves  a  series  of  detailed  steps  that  quickly  overwhelm  the  pencil-and-paper  ap¬ 
proach.  Consequently,  there  is  a  bona  fide  need  for  the  development  of  software  tools 
that  ‘assist’  the  modeler  in  designing  and  implementing  computational  models  of  cel- 
luar  networks.  Two  such  software  packages  are  discussed  in  this  thesis — BioCharon 
and  JigCell. 

In  this  chapter,  1  mainly  focus  on  the  modeling  of  non-biological  and  nat¬ 
urally  occurring  networks.  1  furnish  background  information  on  the  non-biological 
network — the  Brusselator — and  naturally  occurring  networks — the  lysis- lysogeny  path¬ 
way  of  a  mutant  bacteriophage,  the  control  system  that  activates  bioluminescence 
in  Vibrio  hscheri  (V.  fischeri),  molecular  control  of  the  cell  cycle  in  budding  yeast, 
and  a  glycolytic  pathway.  1  also  give  a  brief  description  on  both  software  packages 
to  include  system  configuration,  system  capabilities,  user  interaction,  and  program 
output. 

Besides  describing  each  software  package  and  cellular  network,  1  lay  the  foun¬ 
dation  of  how  computational  models  are  derived  from  cellular  wiring  diagrams.  1 
incorporate  the  rate-reaction  approach  in  each  derivation.  Each  cellular  network 
model  in  this  thesis  is  converted  into  MATLAB  code.  Simulation  results  obtained 
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from  MATLAB  are  regarded  as  baseline  measurements  for  comparison  to  simulation 
results  produced  by  BioCharon  and  JigCell. 

3.2  Rate- Equation  Approach 

As  previously  mentioned  in  section  2.8,  mathematical  modeling  involves  three 
common  approaches:  Boolean,  rate-equation,  and  stochastic.  The  rate-equation 
approach — the  only  approach  that  is  implemented  in  this  thesis — involves  three  main 
steps. 

1.  Convert  a  wiring  diagram  into  a  set  of  reaction  equations. 

2.  Derive  individual  reaction  rate  laws  (i.e.,  reaction  velocities)  from  the  set  of 
reaction  equations  derived  in  step  1. 

3.  Express  each  species  (state  variable)  with  the  appropriate  subset  of  reaction 
rate  velocities  derived  in  step  2. 

Rules  to  derive  individual  reaction  rate  laws  (step  2)  can  be  found  in  [25:8-10]  and 
[12:11-13].  The  Vi’s  in  equations  3.3  and  3.4  are  derived  using  equation  2.7  (just  as 
equation  2.15  results  from  equation  2.7). 

As  a  concrete  example,  I  implement  the  rate-equation  approach  for  the  fol¬ 
lowing  reaction  system  illustrated  in  figure  3.1,  using  steps  one  through  three  and 
basic  laws  of  kinetics  (equations  2.6  and  2.7).  This  reaction  system  is  published  in 
[24:397], 

Step  la: 

ki 

2X  +  At1  3X  (3.1) 

where  k\  is  the  forward  reaction  rate  constant  and  k- \  is  the  reverse  reaction  rate 
constant  for  equation  3.1. 
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Figure  3.1.  This  is  a  possible  wiring  diagram  for  the  system  of  reaction  equations 
proposed  by  [24:397]. 


Step  lb: 


(3.2) 


where  k2  is  the  forward  reaction  rate  constant  and  /c_ 2  is  the  reverse  reaction  rate 
constant  for  equation  3.2. 


Step  2: 


V\  =  — Aq[A"]2[X]  +  k_i[X]3 

(3,3) 

v2  =  —  A^fA]  +  k_2[B] 

(3.4) 

Step  3a: 


d[X] 

dt 


2v\  —  3t'i  +  v2 


=  —V\  +  v2 

=  ki [X]2 [A]  -  [X]3  -  k2 [X]  +  k_2 [B\ 
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Step  3b: 


d[A\ 

dt 


=  v i 


=  —ki[X]2[A\  +  k-i[Xf 


-v2 

k2[X]  -  k_2[B ] 

Upon  inspection,  one  notices  that  ^  =  0.  This  implies  that 

([X}  +  [A]  +  [B])=T1.  (3.5) 

where  Tf  is  a  constant.  Equation  3.5  defines  a  conservation  relation  among  three 
species.  Even  though  concentrations  of  these  species  may  change  over  time,  the  sum 
of  their  concentrations  remain  constant. 

Equation  3.5  can  also  be  rewritten  as  equation  3.6. 

[B]  =  T,  -  (P]  +  [B]).  (3.6) 

Consequently,  the  dynamics  of  this  reaction  system  can  be  expressed  by  a  system  of 
two  non-linear  ODEs  and  one  algebraic  equation.  This  system  is  defined  by  equation 
3.7. 


Step  3c: 


d[B\ 

dt 


«  =  k'lWM  -  V.[A']3  -  k2[X\  +  k_t[B] 

m  =  -kl[x?[A]  +  k-1[x? 

[B]  -  Ti  -  ([A]  +  P]) 
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3.3  Bio  Charon 


BioCharon  is  a  software  package  for  designing,  simulating,  and  analyzing  com¬ 
plex  biomolecular  networks  using  hybrid  systems  (a  hybrid  system  involvs  both  dis¬ 
crete  and  continuous  dynamics)  [5:1].  It  implements  a  graphical  approach  for  user 
interaction  in  building  biomolecular  models.  This  package  is  comprised  of  two  soft¬ 
ware  programs — Bio  Sketch  Pad  and  Charon.  Bio  Sketch  Pad  is  a  GUI-based  drawing 
program  and  Charon  is  a  hybrid-systems  simulation  program  [5:1]. 

The  user  designs  and  simulates  a  biomolecular  model  by  completing  the  fol¬ 
lowing  steps. 

1.  Input  model  into  Bio  Sketch  Pad  (model  definition). 

2.  Translate  model  into  Charon  programming  language. 

3.  Load  translated  model  file  ( modelname.prj )  into  Charon  Visual  Simulator. 

4.  Compile  loaded  model  hie  (select  Generate  Simulator  command). 

5.  Select  simulation  options  ( number  of  steps  and  integration  steps  [i.e.,  integra¬ 
tion  step  size]). 

6.  Select  simulation  type  ( background  simulation  or  display  simulation  trace). 

7.  Start  simulation  (assert  start  command  button). 

Figure  3.2  sketches  out  these  steps. 

3.3.1  Bio  Sketch  Pad.  Bio  Sketch  Pad  (BSP)  is  an  interactive  tool  for 
modeling  and  designing  biomolecular  and  cellular  networks  using  a  graphical  front 
end  [5:2],  This  front  end  (see  figure  3.3)  implements  the  rate-equation  approach  via 
graphical  interaction  with  the  user.  In  this  manner,  the  modeler  defines  a  model  in 
the  GUI  by  assembing  together  defined  nodes  and  arcs.  Additionally,  BSP  automat¬ 
ically  checks  for  syntax  errors  while  the  user  is  building  the  model. 
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Figure  3.2.  This  depicts  the  process  of  how  a  computational  model  is  generated 
from  user  input  via  BioCharon.  The  user  first  designs  a  model  in 
BSP’s  GUI.  Then,  this  completed  model  is  translated  into  the  Charon 
programming  language.  Finally,  this  translated  model  is  loaded  into 
Charon’s  visual  simulator  for  subsequent  simulations. 


Figure  3.3.  This  is  an  entry  into  BSP  that  implements  the  graphical  approach  for 
user  interaction. 
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For  example,  BSP  assures  that  nodes  and  arcs — main  elements  of  BSP — are 
properly  oriented.  Consequently,  the  user  can  only  connect  species  to  reactions  and 
vice  versa.  BSP  cues  the  user  by  coloring  intended  coupling  points.  Green  signifies 
that  the  user  is  attempting  a  legal  operation.  Contrarily,  red  means  that  the  user  is 
attempting  an  illegal  operation. 

The  orientation  and  parameter  definitions  of  nodes  and  arcs  completely  charac- 
erize  the  model.  The  nodes  are  either  species,  chemical  reactions,  or  regulations.  The 
arcs  describe  relations  or  interactions  between  the  nodes  [5:2], 

3.3.2  Charon.  Charon  is  a  hybrid-system  analysis  and  simulation  tool 
[5:1].  After  loading  translated  model  hies  from  BSP  into  the  Charon  visual  simula¬ 
tor,  the  modeler  can  perform  simulations  and  analyze  the  results  within  the  visual 
simulator’s  GUI.  The  process  of  how  a  graphical  model  created  in  BSP  is  translated 
to  Charon  programming  code  is  outlined  in  figure  3.2. 

Charon’s  programming  language  serves  as  an  interface  between  BSP  and  Charon. 
Since  the  modeler  may  not  be  able  to  directly  implement  some  desired  characteristics 
of  a  cellular  network  using  BSP  as  a  front-end  (e.g.,  stochastic  behavior,  inhibition, 
or  Michaelis-Menten  rate  law),  it  may  be  necessary  for  the  modeler  to  include  some 
functions  directly  into  the  Charon  programming  code.  Additional  information  about 
BSP  and  Charon  is  located  in  [5]. 

3.4  JigCell 

JigCcll  is  another  sofware  package  that  can  be  used  for  modeling  cellular  net¬ 
works.  It  is  comprised  of  three  main  software  programs:  JigCcll  Model  Builder 
(JCMB),  JigCcll  Run  Manager,  and  Comparator.  The  user  designs  and  modifies 
models  with  JCMB.  The  Run  Manager  is  the  platform  for  running  simulations  of 
models  created  in  JCMB.  The  user  is  able  to  visually  verify  the  ‘goodness’  of  a  model 
using  the  Comparator.  This  program  generates  a  combined  plot  of  model  predictions 
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and  experimental  temporal  data  points  for  direct  comparison.  Additional  informa¬ 
tion  about  JigCell  is  located  in  [30]. 

3.5  Biomolecular  networks 

“Technological  advancements  combined  with  intensive  DNA  sequencing  ef¬ 
forts  have  generated  an  enormous  database  of  sequence  information  over  the  past 
decade  [21:153].”  Mathematical  modeling  is  used  to  help  bridge  the  gap  between 
this  huge  database  of  sequence  information  and  the  scientific  understanding  of  it. 
Two  types  of  biomolecular  networks  are  generally  modeled  for  the  enhancement  of 
this  understanding — naturally  occurring  and  abstract  networks.  Naturally  occur¬ 
ring  networks  are  specific  biological  systems,  such  as  E.  Coli  and  M.  tuberculosis 
(prokaryotes)  and  Saccharomyces  cervisiae  and  Homo  sapiens  (eukaryotes).  Ab¬ 
stract  networks  are  synthetic  networks.  In  the  remaing  sections  of  this  chapter,  I 
will  give  background  information  on  naturally  occurring  and  abstract  models  as  well 
as  derive  computational  models  from  each  network  wiring  diagram. 

3.6  Non- Biological  Network:  Modified  Brusselator 

“The  accurate  mathematical  description  of  synthetic  networks  provides  the 
foundation  for  describing  complex,  naturally  occuring  networks  [10:268-269].”  In¬ 
sight  gained  from  synthetic  neworks  can  sometimes  be  applied  to  a  class  of  biological 
subsystems.  The  Brusselator  is  a  synthetic  network  that  is  covered  in  this  section. 
“It  is  well  known  that  the  Brusselator  model  reveals  simple  oscillations  and  damped 
oscillations  by  the  suitable  parameters  [18:249].” 

First  proposed  by  Prigogine  and  Lefever,  the  Brusselator  model  is  represented 
by  the  following  system  of  reaction  equations  [22:1697]. 

A  ->  X  (3.7) 
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B  +  X  -»•  Y  +  D 

(3.8) 

2X  +  Y 

-»•  3A 

(3.9) 

X  -»• 

E 

(3.10) 

The  overall  reaction  is  A+B  — *  E+D  [22:1697].  The  concentrations  of  A,  B ,  D,  and 
E  are  assumed  to  remain  constant  and  greater  than  zero.  All  forward  kinetic  rate 
constants  are  set  to  one.  This  system  only  involves  two  intermediate  components 
(species) — X  and  Y  [22:1697].  X  changes  Y  into  X  by  the  autocatalytic  reaction  in 
equation  3.9  [18:249].  The  reaction  rate  for  that  equation  is  increased  as  the  product 
forms.  Furthermore,  the  stoichiometric  coefficients  for  X  in  that  equation  are  two 
and  three  on  the  substrate  side  and  product  side,  respectively. 

Morikawa  et  al.  generalized  equation  3.9.  They  substituted  (2  +  e)  for  the 
stoichiometric  coefficient  for  X  on  the  product  side  of  that  equation.  The  original 
system  of  reaction  equations  is  modified  by  replacing  equation  3.9  with  equation 
3.11. 


2X  +  Y  -»•  (2  +  e)X  (3.11) 

The  free  parameter  e  represents  the  strength  of  the  autocatalytic  process  [18:249]. 
When 

e  is  set  to  one,  equation  3.11  reduces  to  equation  3.9.  I  postulate  that  as  the 
value  of  e  is  increased,  the  steady  state  concentration  of  species  X  is  also  increased. 
Using  steps  two  and  three  in  section  3.2,  the  computational  model  is  derived  from 
the  following  equations. 


Ki  =  -M 

(3.12) 

=  -[Bp] 

(3.13) 
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^3  =  -[A?[F] 

(3.14) 

V4  =  -[X] 

(3.15) 

Reaction  velocities  V\  through  u4  are  collectively  used  to  express  the  dynamics 
(fluxes)  of  species  X  and  Y.  Model  components  A,  B ,  D,  and  E  are  not  regarded  as 
species  or  state  variables  of  the  computational  model,  because  their  concentration 
levels  remain  constant.  The  fluxes  of  species  X  and  Y  are  characterized  by  the 
following  system  of  ODEs. 


d[X] 

dt 


m 

dt 


-Vi  +  v2  +  2v3  -  (2  +  e)v3  +  v4 


-Vi  +  v2  +  2v3  -  2v3  -  ev3  +  u4 


-Vi  +  v2  -  ev3  +  v4 
[A]  -  [B][X]  +  e[Xm  -  [X] 
[A]-([B]  +  l)[X]  +  e[Xm 
-v2  +  v3 

im\  -  [atm 


(3.16) 


(3.17) 


By  setting  [A]  =  a,  [B]  =  b ,  [X]  =  x,  and  [F]  =  y ,  equations  3.16  and  3.17  are 
equivalent  to  the  following  system  of  equations. 

u'T 

—  =  a  —  (b  +  l)x  +  ex2y  (3.18) 

=  bx  —  x2y  (3.19) 

LLL 
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Equations  3.18  and  3.19  can  be  rewritten  as  %  =  /i  and  %  =  /2,  respectively. 
Functions  fi  and  /2  are  defined  by  equations  3.20  and  3.21. 


fi  —  a  —  (b  +  l)x  +  ex2y  (3.20) 

f2  =  bx  —  x2y  (3-21) 

The  critical  point  solution  for  the  system  of  equations  3.16  and  3.17  is  derived 
by  setting  /i  =  0  and  /2  =  0. 

a  —  (6  +  l)x  +  ex2y  =  0  (3.22) 

bx  —  x2y  =  0  (3.23) 


Solving  for  y  in  3.23  results  in  equation  3.24.  Inserting  equation  3.24  into  equation 
3.22  and  then  solving  for  x  results  in  equation  3.25. 


x  = 


6(1  -  e)  +  1 


(3.24) 

(3.25) 


Equations  3.24  and  3.25  coupled  together  imply  y  =  b  ^+1).  Therefore  the 
critical  point  solution  is  represented  by  equation  3.26. 


(x,y) 


6(1  -  e)  +  1 


,b 


'6(1  -e)  +  1' 


(3.26) 


From  equation  3.26,  it  is  clear  that  the  initial  concentration  values  for  species  X  and 
Y  have  no  impact  on  the  critical  point  solution.  It  is  solely  a  function  of  controllable 
parameters  a,  b,  and  e. 
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In  stability  analysis,  time  developments  of  the  fluctuations  (x,  y)  around  the 
critical  point  solution  (ycp  =  ( x,y ))  are  approximated  by  the  following  linearized 
equations  [18:250]. 


dx 

dt 

dy 

dt 


Lux  +  LV2y 


L2\x  +  L22y 


Ln,  L 12,  I/21,  and  L22  are  elements  of  the  Jacobian  matrix  (J).  The  constant  values 
for  these  elements  are  derived  from  the  following  set  of  equations. 


Ln  =  a/liycp)  =  (2e  —  1)6-1 


dx 


t  _  9/i(ycP)  _  ,-2 
Ll2  -  dy  -  ex 


Ln  =  =  -b 

t  —  9/2 (yep)  _  ~2 

^22  -  - Q-y -  -  -X 


These  elements  are  used  to  populate  the  following  Jacobian  matrix. 


1  (2e  -1)6-1 
1  ~b 


At  this  point,  Morikawa  et  al.  ([18:250])  set 


b 


1 

(2e-l)' 


(3.27) 


(3.28) 


Equation  3.28  implies  e  >  1/2,  because  b  >  0.  However,  the  value  of  e  must  be 
an  integer  equal  to  or  greater  than  one  in  order  to  maintain  physical  relevance. 
Subsequent  substitutions  for  b  imply  Ln  =  0  and  L2 1  =  pyvp-  Equation  3.27 
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reduces  to  the  new  Jacobian  matrix  listed  in  equation  3.29. 


0  ex 2 

\  — — - —  — x‘ 2 

V  (2.-1)  x 


(3.29) 


The  trace  [tr(J)],  determinant  [det(J)\,  and  characteristic  equation  [ char(J )] 
are  derived  from  equation  3.29  as  follows. 


tr(J)  = 


L 11  +  T22 

0  H — x2 
'2 


=  —x 


det(J)  = 


L11L22  —  L12L21 
1 


=  0  +  ex 


.(26-1), 


ex 


char(J)  = 


(26  —  1) 

A2  —  tr(J)  A  +  det(J) 


=  A2  +  x2X  + 


ear 


(26  -  1) 


(3.30) 


(3.31) 


(3.32) 


I  use  the  quadratic  formula^  ^  to  solve  for  eigenvalues  A 1  ;2  by  setting 

equation  3.32  to  zero.  In  this  case  a  =  1,  —  b  —  tr(J),  and  c  =  det(J). 


Al,2  — 


tr(J)  ±  \J ( —tr{J ))2  —  4  ( det{J )) 
2 

-£2  ±  ^(£2)2  -  4  (^y) 


-x2  ±  \  x' 


4(<*>) 


(3.33) 


An  observation  from  equation  3.33  is  that  —x2  <  0  (assuming  x  7^  0).  If 
the  expression  a;4  —  4  ( ^e-i) )  <  ^  equation  3.33  is  true,  then  the  values  of 
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A12  contain  negative  real  parts  and  imaginary  parts.  As  a  result,  the  critical  point 
solution  is  stable.  This  expression  can  be  solved  for  x;  however,  it  is  desirable  to 
solve  for  a  controllable  parameter,  such  as  a.  Below,  I  derive  a  condition  that  ensures 
that  the  values  of  Ai^  contain  imaginary  parts. 


a 


x4  —  4 


x1  —  4 


ex 


pe-l) 

e 


.(2c- 1). 

(2e  -  l)2(a)2  4e 


<  0 

<  0 

<  0 


(2e  —  l)2(a)2  <  4e 


2e  —  1 


< 


4e 


2e-  1 J  V(2e-  l)2, 
4e3 


a2  < 


a  <  2 


(2c  -  l)3 


2e  —  1 


(3.34) 


(3.35) 


Equation  3.34  is  derived  by  making  a  substitution  for  x2  (squaring  both  sides  of 
equation  3.25).  The  critical  point  solution  for  the  system  of  equations  3.16  and  3.17 
is  stable  when  conditions  defined  by  equations  3.28  and  3.35  and  e  is  an  integer  equal 
to  or  greater  than  one  are  satisihed. 


3.7  Natural  Networks 

“Although  abstract  models  can  offer  insight  into  basic  mechanisms,  modeling 
must  ultimately  be  connected  to  specific  systems  so  that  verifiable  predictions  can 
be  made  [10:269].”  Since  there  is  a  large  resevoir  of  experimental  data  currently 
available  on  a  diverse  set  of  prokaryotes  and  eukaryotes  [26:247],  the  stage  is  set 
to  enable  tight  coupling  of  candidate  models  to  these  naturally  occurring  systems. 
In  this  section,  three  natural  systems  are  considered  for  modeling  in  JigCell  and 
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BioCharon — bacteriophage,  Vibrio  fischeri,  and  budding  yeast  (Sachharomyces  cere- 
visiae) . 

3.7.1  Lysis-Lysogeny  Pathway.  Hasty  et  al.  propose  a  system  of  reaction 
equations  to  characterize  the  lysis-lysogeny  pathway  of  a  mutant  A  phage.  Unlike 
a  wild  type  A  phage,  the  pathway  of  a  mutant  A  phage  consists  of  two  (instead  of 
three)  operator  sites — OR2  and  OR3.  This  “system  is  a  plasmid  consisting  of  the 
PR  —  Prm  operator  region  and  components  necessary  for  transcription,  translation, 
and  degradation  [11:2076].” 

The  basic  dynamical  properties  of  this  pathway  are  as  follows.  “The  gene 
cl  expresses  repressor  (Cl),  which  in  turn  dimerizes  and  binds  to  the  DNA  as  a 

transcription  factor .  Binding  at  OR2  enhances  transcription,  which  takes  place 

downstream  of  OR3,  whereas  binding  at  OR3  represses  transcription,  effectively 
turning  off  production  [11:2076].”  The  system  of  reaction  equations  that  describe  flux 
of  the  lysis-lysogeny  pathway  of  a  mutant  A  virus  is  described  by  equations  3.36-3.42. 
X,  X2,  and  D  denote  the  concentrations  of  repressor  (a  protein),  repressor  dimer,  and 
DNA  promoter  sites,  respectively.  DX2,  DX2,  and  DX 2X2  denote  binding  to  sites 
OR2,  OR3,  or  both,  respectively.  P  denotes  the  concentration  of  RNA  polymerase, 
and  n  is  the  number  of  proteins  per  mRNA  transcript  [11:2076]. 

Input  Al  (3.36) 

where  r  is  the  basal  rate. 

fci 

2X  kX  X2  (3.37) 

where  k\  is  the  forward  reaction  rate  and  k_\  is  the  reverse  reaction  rate. 

k2 

D  +  X2  DX2  (3.38) 
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where  k2  is  the  forward  reaction  rate  and  k- 2  is  the  reverse  reaction  rate. 

kz 

d  +  x2  C  dx* 

where  ks  is  the  forward  reaction  rate  and  k- 3  is  the  reverse  reaction  rate. 

DX2  +  X2  X  DX2X2 

where  k±  is  the  forward  reaction  rate  and  k- 4  is  the  reverse  reaction  rate. 

DX2  +  P  H  DX2  +  P  +  nX 
where  kt  is  the  reaction  rate. 

x  Ha 

where  k(j  is  the  reaction  rate. 

The  following  reaction  velocities  are  calculated  from  the  previously 
system  of  reaction  equations. 

V\  =  r 

v2  =  -k1[X]2  +  k.1[X2] 
v3  =  -k2[D][X2]  +  k-2[DX2] 
vA  =  -k3[D][X2]+k.3[DX*] 
v5  =  -h[DX 2\[X2]  +  k^[DX2X2\ 
vq  =  -kt[DX2][P] 
v7  =  -kd[X] 


(3.39) 

(3.40) 

(3.41) 

(3.42) 

defined 


3-16 


A  system  of  ODEs — expressing  the  flux  of  each  species — can  now  be  formulated 
from  this  set  of  reaction  velocities. 

=  vi  +  2u2  -  nv6  +  v7 
=  -v2  +  v:i  +  vA  +  v5 

d[D)  . 

dT=V  3  +  ^4 

d\DX  21  ,  , 

=  -V3  +  V5  +  V6  -  Vq 

d\DX*  1 

<it  =  “”4 

dfZ?X2X2l  _  _ 
dt  —  t5 

Since  +  —/ p-  +  =  0  (from  equations  3.45  -  3.48),  there 

exists  a  constant  relation  among  the  concentrations  of  D ,  DX 2,  DX2l  and  DX2X2 
promoter  sites.  At  any  time,  the  total  concentration — denoted  by  dx — remains 
constant  among  these  four  promoter  sites.  The  value  of  dx  is  determined  by  the 
sum  of  initial  concentrations  from  all  four  promoter  sites.  Equation  3.49  defines  that 
constant  relation. 

dT  =  [D]  +  [DX2]  +  [DX*]  +  [DX2X2]  (3.49) 

The  above  system  of  ODEs  reduce  to  the  following  computational  model. 

¥  =  r  +  2(k_1[X2]  -  k^Xf)  +  nkt[DX2][P]  -  kd[X]  (3.50) 

^  =  h[X]2  -  k.,[X2]  -  k2[D][X2]  +  k-2 [DX2\  -  k3[D][X2]+ 

k-:i[DX*2]  -  k4[DX2] [X2]  +  k_A[DX2X2\  (3.51) 

¥  =  k_2[DX2]  -  k2[D][X2]  -  k3[D][X2]  +  k_3[DX*]  (3.52) 

=  k.4[DX2X2\  -  [k.2[DX2]  -  k2[D][X2])  -  k4[DX2][X2]  (3.53) 

=  k3[D][X2]  -  k.3[DX*\  (3.54) 


(3.43) 

(3.44) 

(3.45) 

(3.46) 

(3.47) 

(3.48) 
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d[PX2X2\ 

dt 


k4[DX2][X2]  -  k.4[DX2X2\ 


(3.55) 


Solving  for  [DX2X2]  in  equation  3.49,  a  substitution  is  made  for  [DX2X2]  in 
equations  3.51  and  3.53.  Equation  3.51  reduces  to  equation  3.56,  and  equation  3.53 
reduces  to  equation  3.57. 

^  =  ki[X]2  -  k-i[X2]  -  k2[D][X2]  +  k.2[DX2\  -  k3[D}[X2}  + 

k_z{DX*\  -  k4 [DX2\ [X2]  +  k.4{dT  -([£>]  +  [. DX2 ]  +  [DX*])  (3.56) 

=  k_4(dT  -  ([D]  +  [DX2\  +  [DX*])  -  (k_2[DX2]  -  k2[D][X2])~ 

k4[DX2\[X2\  (3.57) 

As  a  result,  the  computational  model  can  now  be  expressed  by  a  system  of  five  ODEs 
(equations  3.50,  3.56,  3.52,  3.57,  and  3.54).  This  reduced  system  has  a  unique  critical 
point  solution. 

Hasty  et  al.  set  x  =  [A"],  y  =  [X2],  d  —  [D],  u  =  [DX2\,  v  =  [DX g],  and 
z  =  [DX2X2]  [11:2076].  Then,  Hasty  et  al.  reduce  this  system  of  six 
equations  (3.60  -  3.55)  into  two  equations  by  deriving  the  following  system  of  equa¬ 
tions: 


y  =  Kix2 
u  =  K\k2dx2 
v  =  a\K\K2dx2 

z  =  a2(KiK2)2dx4  (3.58) 


where  AT  =  K2  =  ax  =  k3/k2,  and  a2  =  kA/k2. 

Also,  the  total  concentration  of  DNA  promoter  sites  dx  is  constant,  so  that 
dx  —  d  +  u  +  v  +  z.  1  rewrite  this  equation  with  respect  to  d,  so  that  the  entire 
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system  has  only  one  degree  of  freedom. 


kik2x2  +  aikik2(x )2  +  a2{kik2)2(x)A 


(3.59) 


Finally,  Hasty  et  al.  presents  the  nondimensional  equation  that  represents  the 
dynamics  of  the  entire  system. 


/ 

X 


ax2 

1  +  (1  +  ai)x2  +  cr2x 4 


—  'yx  +  1 


(3.60) 


where  a  =  nktpodT/r  and  7  =  kfj/ {r \J~(k\k2). 

The  critical  points  for  3.60  satisfy  the  roots  for  the  fifth  degree  polynomial 
listed  in  equation  3.61.  This  equation  is  obtained  by  setting  the  right  hand  side 
of  equation  3.60  to  zero  and  performing  algebraic  manipulations  on  the  resulting 
equation. 


,y<j2x5  —  a2x4  +  ya;3(l  +  <Ti)  —  x2(a  +  1  +  07)  +  7a:  —  1  =  0  (3.61) 


Five  roots  can  be  determined  from  equation  3.61.  This  is  discussed  further  in  section 
4.4.1. 


3.7.2  Vibrio  fischeri.  As  previously  discussed  in  section  3.4,  V.  fischeri  has 
the  ability  of  bioluminescence.  “The  transcriptional  activation  of  the  lux  genes  in 
the  bacterium  controls  this  luminescence  [1:8]”.  The  set  of  lux  genes  that  controls 
luminescence  is  regarded  as  an  operon. 

An  operon  is  a  sequence  of  closely  associated  genes  that  regulate  enzyme 
production.  An  operon  includes  one  or  more  structural  genes,  which 
carry  information  for  the  synthesis  of  specific  proteins  such  as  enzyme 
molecules,  and  regulatory  sites ,  which  control  the  expresson  of  the  struc¬ 
tural  genes.  A  regulator  gene  works  in  conjunction  with  the  operon  but 
may  be  located  some  distance  from  it  [6:180]. 
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The  lux  regulon  (i.e.,  lux  genes  associated  with  regulation)  is  organized  into  two 
transcriptional  units:  Ol  and  Or.  Ol  contains  the  luxR  gene  and  Or  contains  seven 
genes  (/■moTCDABEG)  [1:8].  These  structural  genes  encode  the  following  proteins. 

luxR  encodes  the  protein  LuxR  [1:8]. 

luxl  “produces  the  protein  LuxI,  which  is  required  for  endogenous  production  of 
the  autoinducer  Ai,  a  small  membrane-permeable  signal  molecule  [1:8].” 

luxA  and  luxB  are  templates  for  luciferase  subunits  [1:8]. 

luxC  and  luxD,  and  luxE  “code  for  proteins  of  the  fatty  acid  reductase,  which 
generates  aldehyde  substrate  for  luciferase  [1:8].” 

luxG  “is  thought  to  encode  a  flavin  reductase  [1:8].” 

The  architecture  of  the  lux  regulon  is  illustrated  in  Figure  3.4.  Additionally, 
that  figure  contains  a  description  of  species  interactions — affecting  luminescence — 
that  takes  place  within  the  cell.  The  pathway  description  is  extracted  from  [1:8]. 

Due  to  the  advanced  nature  of  deriving  the  computational  model  for  charac¬ 
terizing  luminescence  of  V.  fischeri,  I  omit  steps  1-3  as  described  in  section  3.2.  The 
complete  set  of  reaction  equations  is  listed  in  appendix  A  as  input  for  a  JigCell 
calculation.  In  addition,  a  resultant  rate  law  (i.e.,  Uj)  is  automatically  generated  for 
each  reaction  equation  entry  during  model  formulation.  The  following  computational 
model  for  luminescence  is  submitted  by  Alur  et  al.  [1:8].  Descriptions  of  variables 
and  constants  for  this  model  are  listed  in  tables  3.1  and  3.2,  respectively. 


T0  =  kGx o 

Tl  =  Tc(l/;(x 8,  Kco-icdabeg®(cCRP ,  KcRPr,  VcRPr )  +  b)  —  —  kGX\ 

=  Tc($(x 8,  Kco-icdabegi>{cGRPi  KcRPr,  VGRPr )  +  b)  —  H^NA  ~  kGX 2 

x3  =  Ttx  1  -  X3/Hsp  -  rAIRx7x3  +  rCox8  -  KGx3 
x4  =  Tix2-  x4/Hsp  -  kGx 4 
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(3.62) 


x5  =  Tix2  -  x5/Hsp  -  kGx 5 
xti  =  Tix2  -  xA/Hsp  -  kGx6 
x7  =  x0(rAIIx4  -  rAIRx7x3  +  rCox 8)  -  x7/HAI 
x8  =  rAIRx 7x3  -  x8/Hsp  -  rCox8  -  kGx8 

where 

kG  =  kg(  1  -  x0/x0max )  (3.63) 

$(A",  KXm,Vxm)  =  — - -  (3-64) 

K  y  m  +  XVxm 

KXm ,  Vxm)  =  1  —  4>(A",  K Xm i  Vxm)  (3.65) 

is  a  sigmoid  function  that  serves  as  a  switch  for  regulating  (from  low  to  high 
activation  levels)  the  transcription  of  rnRNA.  As  the  concentration  of  the  regulatory 
species  X  increases,  the  value  of  <h  increases.  The  maximum  value  of  $  is  bounded  by 
the  limit  of  <f>,  because  <f>  is  monotone.  In  contrast,  the  ip  function  has  the  opposite 
effect  of  the  sigmoid  function.  As  the  concentration  of  the  regulatory  species  X 
increases,  the  value  of  ip  decreases.  The  minimum  value  is  bounded  by  the  limit  of 
ip  (usually  zero). 


Table  3.1.  Species  Cross  Reference  List  (V.  hscheri)  [1:8 


Species 

Description 

x0 

scaled  population  (population  x  i >b/V) 

X\ 

luxR  (mRNA  transcribed  from  0R) 

X2 

luxICDABEG  (mRNA  transcribed  from  0R) 

x3 

protein  LuxR 

x4 

protein  LuxI 

X5 

protein  LuxA/B  ([LuxA/B]  reflects  cell  luminescence) 

Xq 

protein  LuxC/D/E 

X7 

autoinducer  Ai 

X8 

complex  Co 
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Table  3.2.  Description  of  Model  Parameters  (V.  fischeri)  [1:8 


Parameter 

Description 

Tc 

Maximum  transcription  rate 

Tx 

Maximum  translation  rate 

hrna 

RNA  half-life 

h sp 

Stable  protein  half-life 

Hup 

Unstable  protein  half-life 

Hai 

Ai  half-life 

r  All 

Rate  constant  at  which  LuxI  makes  Ai 

fair 

Rate  constant  Ai  binding  to  LUXR 

f  Co 

Rate  constant  of  Co  dissociation 

VCRPr 

Cooperatively  coefficient  for  CRP 

KcRPr 

Half  maximum  cone,  for  CRP 

Vco—icdabeg 

Cooperatively  coefficient  for  Co 

He  o— iedabeg 

Half  maximum  cone,  for  Co 

b 

Basal  trascription  rate 

vb 

Volume  of  bacterium 

V 

Volume  of  solution 

kg 

Growth  rate 

%0max 

Maximum  population 
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Figure  3.4.  This  diagram — based  on  a  figure  presented  in  [1:8] — depicts  the  reg¬ 
ulatory  control  of  luminescence  in  V.  fishcheri.  Reactions  take  place 
as  follows.  (1)  The  autoinducer  Ai  reversibly  binds  to  protein  LuxR, 
forming  a  complex  Co.  This  complex  subsequently  binds  to  the  lnx 
box.  CRP  (cAMP  receptor  protein)  is  assumed  to  reside  within  the 
cell  at  a  constant  concentration  and  also  binds  to  the  lux  box.  (2) 
The  transcription  (production  of  mRNA  x±)  from  the  InxR  promoter 
(a  regulatory  site)  is  positively  modulated  by  the  binding  of  CRP  to 
its  binding  site.  The  transcription  of  the  InxICDABEG  (production 
of  mRNA  X2)  is  positively  modulated  by  the  binding  of  Co  to  the  lux 
box.  (3)  Growth  in  levels  in  CRP  and  Co  provide  negative  feedback 
(inhibit)  luxICDABEG  and  luxR  transcription,  respectively. 
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3.7.3  Budding  Yeast  Cell  Cycle.  As  stated  in  section  2.1,  a  typical  division 
cycle  of  a  eukaryotic  cell  involves  two  main  phases — S  and  M — temporally  separated 
by  G1  and  G2  phases.  Sometimes,  these  phases  are  distinct  (observable),  and  cell 
division  is  usually  symmetrical.  These  phases  encapsulate  the  following  events.  DNA 
replication  is  initiated,  eventually  forming  sister  chromatids  (S  phase).  These  sister 
chromatids  are  then  aligned  onto  the  metaphase  plate  and  subsequently  segregated 
into  two  nuclei  (M  phase).  Finally,  cell  division  occurs,  resulting  in  one  mother  cell 
and  one  new  daughter  cell.  The  mother  cell  now  resides  in  G1  phase,  pending  the 
initiation  of  another  division  cycle. 

Some  characteristics  of  the  cell  cycle  of  budding  yeast  differ  from  the  behav¬ 
ioral  events  described  previously.  Budding  yeast  cells  divide  asymmetrically  (versus 
symmetrically).  Consequently,  the  mass  of  the  daughter  cell  is  proportionally  smaller 
than  the  mother  cell.  Another  dissimilarity  is  that  budding  yeast  cells  “progress  si¬ 
multaneously  through  S  and  M  phases  (DNA  synthesis,  spindle  formation,  and  chro¬ 
mosome  alignment),  without  any  noticeable  condensation  of  chromosomes  [1:370].” 
Nonetheless,  the  division  cycle  for  eukaryotes  is  closely  regulated  to  ensure  that 
the  correct  order  and  strict  alternation  of  S  and  M  phases  are  enforced  (e.g.,  DNA 
replication  is  complete  before  separation  of  sister  chromatids  is  initiated.). 

“Major  cell  cycle  events  [DNA  synthesis,  bud  emergence,  spindle  formation, 
nuclear  division,  and  cell  separation]  in  budding  yeast  are  controlled  by  a  single 
CDK  [cyclin-dependent  kinase]  (Cdc28)  in  conjunction  with  two  families  of  cyclins: 
Clnl-3  and  Clbl-6  [1:370]”.  CDK  binds  with  one  of  the  cyclin  members  to  form  a 
heterodimer  (e.g.,  Clb5/Cdc28).  This  binding  ‘activates’  the  complex  and  enables 
it  to  carry  out  its  regulatory  activities.  Chen  et  al.  list  the  functions  of  each  het¬ 
erodimer  [8:370]. 

•  Clnl/Cdc28  and  Cln2/Cdc28  influence  budding  and  spindle  pole  body  dupli¬ 
cation. 
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•  Cln3/Cdc23  influences  the  size  at  which  newborn  cells  execute  Start. 

•  Clb5/Cdc28  and  Clb6/Cdc28  influence  timely  DNA  replication. 

•  Clb3/Cdc28  and  Clb4/Cdc28  influence  DNA  replication  and  spindle  formation. 

•  Clbl/Cdc28  and  Clb2/Cdc28  influence  proper  completion  of  mitosis. 

Novak  partitions  the  above  heterodimers  into  two  different  sets.  CycB/Cdc28  is  the 
set  that  contains  a  Clb  part  [20:2],  Cln/Cdc28  is  the  set  of  that  contains  a  Cln 
part  [20:4], 

Cell  division  events  are  regulated  by  the  interaction  of  CycB/Cdc28  with 
what  Novak  terms  as  enemies  and  helper  molecules.  There  are  two  types  of  helper 
molecules — helpers  of  the  enemies  and  helpers  of  CycB/Cdc28.  The  enemies  and 
their  helper  molecules  collectively  inhibit  the  activity  of  CycB/Cdc28.  On  the  other 
hand,  helper  molecules  can  reinitiate  Cdk/CycB  activity  (when  a  certain  condition  is 
satisfied),  and  both  eventually  inactivate  enemies  and  their  friends.  Table  3.3  sum¬ 
marizes  various  molecular  interactions  that  modulate  CycB/Cdc28  activity  [20:2-4], 


Table  3.3.  Cdk/CycB  Molecular  Interactions 


Name 

Relation  Type 

Interaction  Effect 

Cdhl 

Enemy 

Promotes  Cdk/CycB  degradation 

Sicl 

Enemy 

Binds  to  Cdk/CycB,  causing  inactivation 

Cdc20 

Enemy  helper 

Promotes  Cdk/CycB  degradation 

Cdk/Cln  kinases 

Cdk/CycB  helper 

Assists  in  initiating  Start  transition 

Due  to  the  antagonisic  interaction  between  CycB/Cdc28  and  enemies,  the  con¬ 
trol  system  for  cell  division  operates  from  two  alternative  states — State i  and  State 2. 
CDK  activity  (combined  activities  of  Gib’s  and  Gin’s)  is  low  in  State i,  while  their 
enemies — Sicl  and  Cdhl  are  high.  Therefore,  concentration  levels  of  the  Gib’s  and 
Gin’s  remain  mostly  low  during  State i.  On  the  other  hand,  CDK  activity  is  high  in 
State 2.  As  a  result,  the  concentrations  of  the  Gib’s  and  Gin’s  are  high,  where  their 
enemies’  concentrations  remain  mostly  low  during  State 2.  Along  these  lines,  Chen 
et  al.  cite  the  following  quote  [8:370]. 
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Nasmyth  (1996)  has  proposed  that  the  heart  of  the  budding  yeast  cell 
cycle  is  an  alternation  between  two  self-maintaining  states:  the  G1  state, 
in  which  APC  is  active,  CDK  activity  is  low,  and  orgins  are  licensed; 
and  the  S/M  state,  in  which  APC  is  shut  off,  CDK  activity  is  high,  and 
origins  are  fired  and  incapable  of  bring  again. 

Control  is  transferred  between  self-maintaining  states  by  means  of  enabled 
transitions.  The  budding  yeast  cell  cycle  has  two  major  transitions — Start  and 
Finish.  Start  is  enabled  when  the  cell  reaches  critical  mass.  Once  Start  is  en¬ 
abled,  Cln/Cdk  helper  molecules  initiate  the  buildup  of  CycB/Cdc28  molecules  by 
activating  transcription  factors  that  synthesize  CycB  molecules  (the  concentration 
of  Cdc28  is  assumed  to  exist  in  excess).  The  combined  concentration  of  CycB/Cdc28 
and  Cln/Cdc28  molecules  eventually  increase  to  a  level  that  allows  them  to  inibit  the 
activity  of  the  enemies  and  their  helper  molecules.  At  this  stage,  the  environment  is 
primed  for  a  new  round  of  DNA  replication. 

Finish  is  enabled  when  DNA  is  fully  replicated  and  chromosomes  are  aligned 
on  the  metaphase  plate.  Once  Finish  is  enabled,  Cdc20 — helper  molecule  for  the 
enemies — is  activated.  Indirectly  Cdc20  promotes  dissociation  of  sister  chromatids, 
activation  of  Hctl  (partly  responsible  for  the  degradation  of  Clb2/Cdc28),  and  ac¬ 
tivation  of  Swi5  (a  transcription  factor  that  initiates  the  synthesis  of  Sicl)  [8:373]. 
Ultimately,  increasing  activity  levels  of  the  enemies  and  their  helper  molecules  nul¬ 
lify  CycB/Cdc28  activity.  Figure  yeastcyclc  is  the  molecular  model  for  control  of 
CDK  activities.  It  is  extracted  from  the  consensus  model  published  by  Chen  et  al. 
[8:373].  Chen  et  al.  list  the  following  events  that  occur  when  either  transition  is 
initiated.  Detailed  information  on  the  cell  division  cycle  for  budding  yeast  can  be 
found  in  [8:369-391], 

Due  to  the  advanced  nature  of  deriving  the  computational  model  for  charac¬ 
terizing  system  control  for  the  cell  division  cycle  of  budding  yeast,  I  omit  steps  1-3  as 
described  in  section  3.2.  The  complete  set  of  reaction  equations  is  listed  in  appendix 
A.  Chen  et  al.  present  the  following  system  of  ODEs  that  characterizes  cell  division 
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Figure  3.5.  CDK  heterodimers  are  represented  by  their  cyclin  partners  (e.g., 
Clb3/Cdc28  =  “Clb3”),  and  redundant  cyclins  are  listed  as  {Clnl  +  Cln2 
=  “Cln2”}  and  {Clb5  +  Clb6  =  “Clb5”}.  At  Gl,  the  cell  has  few  cyclin 
molecules.  When  the  cell  reaches  critical  mass,  Cln3/Cdc28  and  Bck2  acti¬ 
vate  transcription  factors  (TF)  SBF  and  MBF  by  phosphorylation  [8:373]. 
Cln2  and  Clb5  begin  to  accumulate.  Clb5  accumulates  as  inactive  trimers 
of  Clb5/Cdc28/Sicl  [8:373].  Rising  levels  of  Cln2/Cdc28  activity  initiate 
budding,  inactivate  Hctl,  and  assist  in  phosphorylating  Sicl  for  degrada¬ 
tion.  “When  Sicl  is  destroyed,  Clb5/Cdc28  activity  rises  abruptly  and 
drives  the  cell  into  S  phase  [8:373].”  Since  Sicl  is  destroyed,  Clb2/Cdc28 
is  able  to  activate  its  own  TF  (Mcrnl).  Rising  levels  of  Clb2/Cdc28  pro¬ 
vide  negative  feedback  to  SBF.  As  a  result,  Clb5/Cdc28  level  starts  to  fall. 
Clb2/Cdc28  induces  the  cell  to  progress  through  mitosis  (M  Phase).  Cdc20 
is  activated  when  chromosomes  are  aligned  and  DNA  synthesis  is  complete. 
Sister  chromatids  are  then  separated,  Hctl — that  promotes  Clb2/Cdc28 
degradation — and  Swi5  are  activated.  Clb’s  and  Cln’s  are  eventually  de¬ 
stroyed  or  inactived.  Sicl  levels  continue  to  rise  via  an  active  Swi5,  and  the 
cell  returns  to  Gl  [8:373]. 
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control  for  budding  yeast  [8:374], 


ft[ Cln2 }  =  (k'sn2  +  k”n2 [SBF})mass  -  kd,n2[Cln2 ] 

^[Clb2]T  =  {k's  b2  +  k"sb2 [M cml])mass  -  Vd)b2[Clb2}T 
£t[Clb5\T  =  {k'sb 5  +  k”sbb[MBF])mass  -  Vdtb5[Clb5\T 

ftisic  1]t  =  k's,c  1  +  *a,cl[^*5]  -  (fcdl,cl  +  Jd2  J+gcl]T  [^cl]T 

|[C/62/5zcl]  =  kas,b2[Clb2][Sicl\  -  (kdi,b2  +  Vd,b2  +  Kdl,cl  +  jd2^[Sici]T)  [Clb2/Sic  1] 
ft[Clb5/Sicl]  =  kas,b2[Clb5}[Sicl}  -  (kdi,b s  +  +  7Qi,ci  +  jd2J+isici]r)  [C^/Sicl] 

|[Cdc20]T  =  (k'St20  +  k">20[Clb2])  -  kd,2O[Cdc20\T 
^[Cdc20\  =  ka,2O([Cdc20\T  -  [Cdc20])  -  (Vit 20  +  kdj20)[C  dc20\ 

d.  rc7„+il  _  (K  tl+k”M[Cdc20])(Hctl]T-[Hctl])  ViM[ Hctl] 
dt[nciL\  -  JaM+[Hctl\T-[Hctl\  Ji:tl+[Hctl\ 

■^[mass]  =  n{mass) 

2 ]j[ORI ]  =  kS}0ri([Clb5\  +  £ori,b2[C'lb2])  —  kd^ori[ORI] 

[BU D]  =  kS)bud([Cln2]  +  [C7n3]*  +  £bud,bh[Clbh])  -  kd,bud[BUD\ 

(3.66) 


where 


[Clb2}T  =  [Clb2]  +  [Clb2/Sicl] 

[Clbh]T  =  [Clb5]  +  [Clbb/Sicl] 

[5icl]r  =  [Szcl]  +  [Clb2/ Sicl]  +  [C765/5icl]  (3.67) 


3.7.4  Glycolysis.  “Glycolysis  is  the  metabolic  pathway  used  by  most  au¬ 
totrophic  and  heterotrophic  organisms,  both  aerobes  and  anaerobes,  to  begin  to 
break  down  glucose  [6:114].”  In  the  Embden- Meyerhof  pathway,  the  breakdown  of 
glucose  occurs  in  10  steps  or  catalyzed  reactions  [6:A27].  These  10  steps  encompass 
three  stages: 
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1.  Confine  substrate  (glucose)  inside  the  cell  and  form  phosphorylated  six-carbon 
units  (fructose  1,  6-diphosphate). 

2.  A  six-carbon  unit  is  split  into  two  three-carbon  units,  and  two  ATP  molecules 
are  produced.  Two  NAD  molecules  are  reduced. 

3.  Pyruvic  acid  is  formed,  and  two  more  ATP  molecules  are  produced.  [6:A27] 

Figure  A  illustrates  this  metabolic  pathway. 

Each  reaction  (step)  in  a  metabolic  pathway  is  controlled  by  a  particular  en¬ 
zyme.  Ten  enzymes  in  the  Embden-Meyerhof  pathway  are  the  terms  in  figure  A 
that  contain  the  suffix  ‘ase’.  As  previously  stated  in  section  2.3.4,  a  single  catalyzed 
reaction  can  be  divided  into  two  separate  reactions.  In  that  section,  equations  2.10 
and  2.11  are  derived  from  equation  2.9.  Equations  2.13  and  2.14  are  derived  from 
equation  2.12.  Along  these  lines,  the  set  of  10  catalyzed  reactions  is  represented 
by  a  set  of  20  equivalent  reaction  equations.  Table  3.5  lists  the  set  of  10  reaction 
equations  derived  from  figure  A.  Table  3.4  indicates  the  correspondence  between 
the  entries  depicted  in  figure  A  and  the  |/j’s  (for  i  —  1  to  37)  which  appear  in  ta¬ 
bles  3.5,  3.6,  3.7,  and  equation  3.68.  Table  3.6  lists  the  resultant  set  of  20  reaction 
equations  that  is  derived  from  table  3.5.  Each  pair  of  reaction  equations  in  table  3.6 
corresponds  to  a  single  reaction  in  table  3.5.  For  example,  reactions  1A  and  IB  in 
table  3.6  correspond  to  reaction  1  in  table  3.5,  and  reactions  2A  and  2B  correspond 
to  reaction  2,  and  so  on. 

Each  CatTrarij  entry  (for  j  —  1  to  10)  in  table  3.4  is  an  enzyme-substrate 
complex  that  corresponds  to  a  reaction  in  table  3.5.  For  each  one-way  reaction  in 
table  3.5,  each  CatTrarij  is  analogous  to  ES  in  equations  2.10  and  2.11.  For  each 
two-way  reaction  in  table  3.5,  each  CatTrarij  is  analogous  to  ES  in  equations  2.13 
and  2.14.  The  ty.’s  (for  k  =  1  to  20),  which  appear  in  equation  3.68,  are  defined  in 
table  3.7. 
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Alpta-D-Gkose 


Mg++,  Hexokinase 


Alplii-D-Gkose-6-Phosptate 


Figure  3.6.  The  glycolytic  pathway  involves  10  catalyzed  reactions.  As  a  result 
glucose  is  broken  down  to  pyruvate  [6:A27]. 
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I  derive  the  following  system  of  ODEs  from  tables  3.4,  3.6,  and  3.7. 


d\yi] 

dt 

d\V2] 

dt 

t%3] 

dt 

d[y±\ 

dt 

d\ys] 

dt 

d[ye) 

dt 

<%7] 

dt 

d[ys\ 

dt 

djya} 

dt 

d[yio\ 

dt 

d[yu\ 

dt 

d[yi2] 

dt 

^ryi3l 

dt 


=  V4 

=  Vi  +  V5  -  V14  -  V20 
=  ~V2  -  Vq  +  V 13  +  Vi9 

=  Vi  -  V2  +  v5  -  V6  +  Vi3  -  V\4  +  Vn  -  Vis  +  v19  -  v20 
=  Vi  -  v> 

=  —Vx  +  v2 
=  -V2  +  n3 
=  V3  -  v4 
=  -V3  +  n4 
=  —v4  +  n5 
=  V5  -  v6 
=  -v5  +  n6 
=  -VQ  +  v7 


d[yi4,] 

dt 

d[y  w] 
dt 

d[yi6] 

dt 

d\yi7 1 
dt 

d[yis] 

dt 

d[yi9] 

dt 

d\y  2o\ 
dt 

d\y2i] 

dt 

d[y22] 

dt 


=  V7  -  Vs 

=  -v7  +  v8 

=  -Vs  -  Vio  +  nil 

=  —vs  +  n9 

=  Vg  -  VW 
=  -Vg  +  Vio 
=  VU 

=  nil 
=  -ni2 
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d[y  23] 

dt 

d[y 24] 
dt 

d\y 25] 
dt 

d\y2o\ 

dt 

d[y2i] 

dt 

d[y28] 

dt 

d\y29] 

dt 

d\y  30] 
dt 

d[y3i] 

dt 

d[y  32] 
dt 

d\y33] 

dt 

d\y34] 

dt 

d[y  35] 
dt 

d[y  36] 
dt 

d\y3i] 

dt 


"11  ~  "12 
-"11  +  v12 
~Vl2  +  "13 
"13  -  "14 
“"13  + 

— "14  +  "l5 
"l5  ~  '-’Hi 
"15  +  "16 
-"16  +  "17 
"17  -  "18 

—"17  +  Vis 

-"18  +  "19 
"19  —  "20 

-"19  +  V20 
-"20 


(3.68) 


3-32 


Tabic  3.4.  Species  Cross  Reference  List  (Glycolytic  Pathway) 


Variable  Species 


Alpha-D-Glucose 


ATP 


ADP 


Mg++ 


Hexokinase 


CatTrani 


Alpha-D-Glucose- 6-phosphate 


Phosphoglucoisomerase 


CatTran2 


Beta-D-Fructose-6-phosphate 


Phosphofructokinase 


CatTrang 


Beta-D-Fructose-1,  6-bisphosphate 


Aldolase 


CatTrani 


Glyceraldehyde-  3-phosphate 


Dihydroxyacetone  phosphate 


Triose  phosphate  isomerase 


CatTrang 


NAD+ 


Pi 


NADH 


Glyceraldehyde-3-phosphate  dehydrogenase 


CatTrang 


1,  3- Bisphosphogly cerate 


Phosphoglycerate  kinase 


CatTrani 


3-Phosphoglycerate 


Phosphoglycerate  mutase 


CatTrang 


2-Phosphoglycerate 


Enolase 


CatTrang 


Phosphoenolpyruvate 


Pyruvate  kinase 


CatTranw 


Pyruvate 
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Table  3.5.  Reaction  Equations  (Glycolytic  Pathway) 


Reaction 

Reaction  Equation 

1 

2/l  +  Vi  +  2/4  +  2/5  -►  2/3  +  2/4  +  2/5  +  2/7 

2 

V"  +  Vs  ^  2/8  +  l/io 

3 

2/2  +  2/4  +  2/10  +  2/ii  — ^  2/3  +  2/4  +  2/ii  +  V l 3 

4 

2/13  +  2/14  ^  2/14  +  2/16  +  2/17 

5 

1/17  +  2/18  ^  2/16  +  2/18 

6 

2/16  +  1/20  +  1/21  +  1/23  ^  1/22  +  1/23  +  1/25 

7 

1/3  +  1/4  +  1/25  +  1/26  ^  1/2  +  1/4  +  1/26  +  1/28 

8 

1/28  +  1/29  ^  1/29  +  1/31  +  V\1 

9 

1/4  +  1/31  +  1/32  +  1/23  ^  2/4  +  1/32  +  1/34 

10 

2/3  +  2/4  +  1/34  +  1/35  — ^  2/2  +  1/4  +  1/35  +  1/37 
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Table  3.6. 


Reaction  Equations  with  Catalytic  Transitions  (Glycolytic  Pathway) 


Reaction 

Reaction  Equation 

1A 

ki 

2/l  +  2/2  +  2/4  +  2/5  *2  2/6 

IB 

2/6  ^  2/3  +  2/4  +  2/5  +  2/7 

2A 

2/7  +  2/8  *5  2/9 

2B 

*6 

2/9  ^  2/8  +  2/io 

3A 

kg 

2/2  +  2/4  +  2/io  +  2/n  *9  2/i2 

3B 

2/12  ^  2/3  +  2/4  +  2/11  +  2/13 

4A 

*12 

2/13  +  2/14  *T3  2/15 

4B 

k  14 

2/15  *16  2/14  +  2/16  +  2/17 

5A 

*16 

2/17  +  2/18  *17  2/19 

5B 

*18 

2/19  *19  2/16  +  2/18 

6A 

*20 

2/16  +  2/20  +  2/21  +  2/23  *21  2/24 

6B 

*22 

2/24  *23  2/22  +  2/23  +  2/25 

7A 

*24 

2/3  +  2/4  +  2/25  +  2/26  *25  2/27 

7B 

*26 

2/27  *27  2/2  +  2/4  +  2/26  +  2/28 

8A 

*28 

2/28  +  2/29  *29  2/30 

8B 

*30 

2/30  *31  2/29  +  2/31 

9A 

*32 

2/4  +  2/31  +  2/32  *33  2/33 

9B 

*34 

2/33  *35  2/4  +  2/32  +  2/34 

10A 

*36 

2/3  +  2/4  +  2/34  +  2/35  *37  2/36 

10B 

2/36  ^  2/2  +  2/4  +  2/35  +  2/37 
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Table  3.7. 


Derived  Reaction  Velocities  for  Emb den- Meyerhof  Glycolytic  Pathway 


Reaction 

Velocity 

1A 

vi  =  -kiyiy2y^  +  k2y6 

IB 

v2  =  -k3y6 

2A 

v3  =  -k±y7y8  +  k5y9 

2B 

va  =  - hy9  +  k7y8y  w 

3A 

=  — ^82/22/42/io2/n  +  k9yi2 

3B 

vq  =  —  kioyi2 

4A 

V7  —  —knyi3yi4  +  ki2yi§ 

4B 

v8  =  —k\3yi5  +  kuyuyioyn 

5A 

v9  =  —  ki5ynyis  +  ki6yi9 

5B 

V10  =  — ^172/19  +  ^182/162/18 

6A 

Vu  =  —ki9yioy29y2iy23  +  k29y24 

6B 

V{2  =  —k2\y24  +  k22y22y23y23 

7A 

Vl3  =  —  ^232/32/42/252/26  +  ^24?/27 

7B 

^14  =  —k23y27  +  ^262/22/42/262/28 

8A 

^15  —  —k27y28y29  +  k28y3o 

8B 

?’i6  =  ~k29y3o  +  k30y29y3i 

9A 

V17  =  —  k3iy4y3iy32  +  k32y33 

9B 

vis  =  — k33y33  +  k34y\y32y34 

10A 

v  1 9  =  —k35y3y4y34y35  +  k36y3e 

10B 

V20  =  k37y3o 

3-36 


IV.  Computational  Model  Results 

4-1  Overview 

In  chapter  III,  I  discuss  the  following  reaction  models:  modified  Brusselator, 
lysis-lysogeny  pathway  of  a  mutant  bacteriophage,  quorum  sensing  system  that  con¬ 
trols  luminescence  in  V.  fischeri,  control  system  for  cyclin  dependent  kinases  (CDKs) 
in  the  budding  yeast  cell  cycle,  and  a  glycolytic  pathway.  One  way  of  shedding  light 
on  the  overall  functionality  of  these  reaction  systems  is  to  derive  a  computational 
model  (system  of  ODEs)  for  each  reaction  system.  Computational  models  are  gen¬ 
erally  used  to  track  how  state  variables  (concentrations  of  each  species)  evolve  over 
a  specified  time  interval.  The  rate-equation  approach  for  deriving  computational 
models  is  the  method  exclusively  used  in  this  thesis.  This  approach  involves  the 
following  steps. 

1.  Convert  a  wiring  diagram  into  a  set  of  reaction  equations  (e.g.,  table  3.6). 

2.  Derive  individual  reaction  rate  velocities  from  the  set  of  reaction  equations 
derived  in  step  1  (e.g.,  table  3.7). 

3.  Express  each  state  variable  (species)  with  the  appropriate  subset  of  reaction 
rate  velocities  derived  in  step  2  (e.g.,  equation  3.68). 

BioCharon  and  JigCell  are  software  tools  that  serve  as  aids  in  building  biomolec- 
ular  models,  deriving  (via  rate-equation  approach)  computational  models,  and  simu¬ 
lating  these  models.  Steps  1-3  are  automated  in  BioCharon.  Steps  2-3  are  automated 

in  JigCell. 

In  this  chapter,  I  present  general  simulation  results  that  highlight  certain 
qualitative  features  (e.g.,  conservation  relation  among  species)  of  the  previously 
mentioned  reaction  models.  I  then  compare  simulation  results  from  MATLAB, 
BioCharon,  and/or  JigCell.  Metrics  for  these  comparisons  are  explained  in 
section  4.2. 


4-1 


Simulation  results  are  obtained  by  MATLAB,  BioCharon,  and  JigCcll  for  the 
modified  Brusselator  model  and  lysis-lysogeny  pathway.  Simulation  results  for  the 
quorum  sensing  system  are  obtained  by  MATLAB  and  JigCell.  This  system  is  not 
simulated  using  BioCharon,  because  the  cost  is  extremely  high  in  simulation  time 
and  memory  (an  extremely  small  integration  step  size  is  required  to  obtain  accurate 
results).  Simulation  results  for  the  control  system  of  CDKs  in  the  budding  yeast  cell 
division  cycle  are  obtained  from  MATLAB  and  JigCcll.  An  instance  of  this  model  is 
not  created  using  BioCharon,  because  the  user  is  not  able  to  explicitly  define — via 
Bio  Sketch  Pad — key  aspects  of  this  model.  For  example,  the  user  is  not  able  to 
explicitly  define  Michaelis-Menten  rate  laws  and  functions  for  transcription  factors 
(e.g.,  Vd,b 2  and  in  equation  3.66).  Simulation  results  for  the  glycolytic  pathway 
are  soley  obtained  by  MATLAB,  because  the  computational  time  for  MATLAB  is 
substantially  less  than  that  for  both  BioCharon  and  JigCell. 

4-2  Experimental  Design 

All  or  some  of  the  following  metrics  are  used  in  comparing  simulation  results: 
absolute  error,  sum  of  absolute  errors,  and  relative  error.  Absolute  error  is  the 
absolute  value  of  the  difference  between  simulation  results  obtained  from  MATLAB 
and  either  BioCharon  or  JigCell.  This  metric  is  used  in  all  comparisons. 

Sum  of  absolute  errors  is  the  sum  of  individual  absolute  errors  among  all  state 
variables  (species)  for  a  specified  trial.  This  metric  is  used  during  the  comparison 
of  simulation  results  for  the  lysis-lysogeny  pathway  in  order  to  ascertain  overall 
accuracy  for  this  reaction  system.  This  metric  is  used  soley  in  table  4.9. 

Relative  error  is  used  in  comparing  simulation  results  for  the  quorum  sensing 
system.  This  additional  metric  is  used  because  simulation  results  among  species 
vary  several  orders  of  magnitude.  Absolute  error  among  species  having  significantly 
larger  values  is  more  likely  to  have  larger  absolute  error  than  among  species  with 
significantly  smaller  values.  For  example,  luxICDABEG  has  a  value  of  9.80402,  but 
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LuxR  has  a  value  of  5208.60414  in  table  4.12.  Relative  error  is  used  to  better  gauge 
the  degree  of  error  among  all  species.  Relative  error  is  ^MATL\M^rLABRth^t r^Re3'ult^  ■  This 
metric  is  used  in  tables  4.11,  4.12,  4.13,  4.14,  4.15,  and  4.16. 

“With  a  user  community  more  than  500,000  strong  spread  throughout  indus¬ 
try,  government,  and  academia,  MATLAB  is  the  recognized  standard  worldwide  for 
technical  computing  [16:2]”.  In  this  regard,  simulation  results  obtained  from  MAT- 
LAB  are  regarded  as  baseline  results  (i.e.,  error  =  0).  Therefore,  the  accuracy  of 
simulation  results  generated  from  BioCharon  and  JigCcll  are  gauged  from  simulation 
results  generated  by  MATLAB. 

System  parameters  for  BioCharon,  JigCell,  and  MATLAB  are  configured  in 
the  following  manner.  The  integration  step  size  for  the  ordinary  differential  equation 
(ODE)  solver  of  BioCharon  is  set  to  0.01  during  all  simulations.  I  performed  several 
pilot  studies  in  order  to  decide  which  step  size  to  use.  A  smaller  step  size  may 
slightly  increase  accuracy  (depending  on  the  system);  however,  total  simulation  time 
is  greatly  increased  as  the  integration  step  size  is  shortened.  The  ODE  solver  for 
JigCell  is  set  to  ‘stiff’  during  all  simulations  (a  solver  that  is  specifically  designed  to 
approximate  solutions  for  a  stiff  system  of  differential  equations).  The  ODE  solver 
for  MATLAB  is  set  ‘ode23s’  at  the  default  accuracy  (relative  tolerance  of  IE-3  and 
vector  of  absolute  error  tolerances  of  IE-6  by  default).  The  system  parameter  that 
adjusts  total  simulation  time  is  synchronized  for  all  software  tools. 

In  obtaining  simulation  results,  each  of  the  five  computational  models  is  sub¬ 
mitted  as  a  workload  to  the  appropriate  software  tool.  Each  workload  submission 
on  either  MATLAB,  BioCharon,  and  JigCell  has  an  accompanying  set  of  initial  con¬ 
ditions  (species  concentrations  and  model  parameters).  Initial  concentrations  for 
species  are  varied  in  sections  4.3,  4.4,  4.5,  and  4.6.  While  maintaining  the  same 
initial  concentrations  for  species,  a  single  model  parameter  is  varied  in  section  4.7. 
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4-3  Modified  Brusselator  Model 

The  Brusselator  model  “reveals  simple  oscillations  and  damped  oscillations  by 
the  suitable  parameters  [18:249].'"  I  describe  a  modified  version  of  the  Brusselator 
model  in  section  3.6.  The  modification  involves  the  addition  of  a  free  parameter 
e.  This  parameter  reveals  the  strength  of  the  autocatalytic  process  for  species  X 
(equation  3.11)  [18:249].  If  e  is  assigned  an  integer  value  greater  than  or  equal 

to  one  and  model  parameters  a  and  b  are  defined  using  equations  3.35  and  3.28, 
respectively,  the  critical  point  (x,y) — defined  by  equation  3.26 — is  stable. 

Figure  4.1  illustrates  10  critical  points  for  the  modified  Brusselator  system. 
Each  critical  point  is  obtained  by  using  equation  3.26  and  imposing  the  above  con¬ 
straints  for  parameters  a,  b  and  e.  The  integer  value  of  e  is  varied  from  1  to  10.  It  is 
clear  from  figure  4.1  that  the  critical  point  concentration  for  species  x  increases  as  e 
is  increased.  This  behavior  is  consistent,  because  e  serves  as  the  autocatalytic  power 
control  for  x.  It  follows  that  as  the  level  of  e  is  increased,  the  operating  level  of  x  is 
increased.  Figure  4.1  also  illustrates  that  the  relative  effect  of  e  on  x  diminishes  as 
e  is  increased. 

4-3.1  General  Results.  Figures  4.2  and  4.3  illustrate  numerical  simulations 
for  the  Brusselator  system  performed  on  MATLAB.  The  value  for  parameters  a  and 
b  are  selected  in  accordance  with  the  previously  mentioned  constraints  that  ensure 
a  stable  critical  point.  In  both  cases,  a  =  0.7,  b  =  0.5,  and  e  =  1.  In  figure  4.2,  the 
initial  concentrations  of  x  and  y  are  zero.  In  figure  4.3  the  initial  concentrations  for 
both  species  are  five.  Despite  the  significant  difference  between  initial  conditions, 
the  resulting  trajectory  is  back  to  the  same  critical  point:  x  =  0.7  and  y  ~  1.4. 
This  behavior  of  both  trajectories  returning  to  the  same  critical  point  is  illustrated 
in  figure  4.4  (phase  plane  plots).  The  behavior  of  returning  back  to  a  critical  point 
after  applying  large  perturbations  away  from  that  critical  point  alludes  to  the  possible 
existence  of  a  global  critical  point  for  the  modified  Brusselator  reaction  system. 
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Critical  Point  Concentration 


Figure  4.1.  As  e  is  increased  (with  a  =  0.7),  the  critical  point  solution  for  species 
x  increases  and  y  decreases.  It  appears  that  as  e  approaches  10,  the 
critical  point  solution  approaches  x  =  1.33  and  y  =  0.04. 

4-3.2  Comparison  Results.  A  set  of  three  simulations  are  performed  using 
MATLAB,  BioCharon,  and  JigCell.  The  parameters  are  a,  b ,  x,  and  y.  and  the 
single  factor — a  parameter  that  is  varied  in  each  trial — is  e.  The  values  of  a  and  b 
are  selected  to  ensure  a  stable  critical  point.  For  all  three  trials,  the  concentration 
of  a  is  held  constant  at  0.7,  and  the  initial  concentrations  of  x  and  y  are  zero.  The 
constant  concentration  of  b  is  dependent  on  the  single  factor  e  (via  equation  3.28). 
This  factor  is  assigned  integer  values  1  through  3.  These  parameter  and  factor  values 
are  listed  in  table  4.1. 

Table  4.2  lists  simulation  results  from  MATLAB,  BioCharon,  and  JigCell  for 
species  x  and  y.  These  results  are  critical  point  approximations  (i.e.,  the  rate  of 
change  for  all  state  variables  is  zero)  for  the  modified  Brusselator  system.  The 
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4.2.  Parameter  settings  for  the  modified  Brusselator  model  are 
b  =  0.5,  and  e  —  1.  Initial  concentrations  are  x  =  0  and  y  = 


concentration 


0  10  20  30  40  50  60  70  80  90  100  110 


[time  units] 


[time  units] 


Figure  4.3.  Parameter  settings  for  the  modified  Brusselator  model  are  a 
b  =  0.5,  and  e  =  1.  Initial  concentrations  are  x  =  5  and  y  —  5. 
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Figure  4.4.  Parameter  settings  for  the  modified  Brusselator  model  are  a  =  0.7, 
b  =  0.5,  and  e  —  1.  The  initial  concentrations  for  the  top  plot  are 
x  =  0  and  y  —  0,  and  the  initial  concentrations  for  the  bottom  plot 
are  x  =  5  and  y  —  5. 


4-8 


results  in  that  table  are  separated  by  trial  numbers  1  through  3.  Critical  point 
approximations  for  each  software  tool  are  identical:  x  =  1.42857  and  y  =  1.42857 
for  trials  1  through  3.  Therefore,  absolute  error  for  all  critical  point  approximations 
is  zero. 


Table  4.1.  Initial  Conditions  (Modified  Brusselator  Model) 


Tidal 

A 

B 

X 

Y 

e 

1 

0.7 

1.0 

0 

0 

1 

2 

0.7 

l 

3 

0 

0 

2 

3 

0.7 

0.2 

0 

0 

3 

Table  4.2.  Critical  Point  Results  (Modified  Brusselator  Model) 


Tidal 

Application 

X 

Y 

1 

MATLAB 

0.70000 

1.42857 

BioCharon 

0.70000 

1.42857 

JigCcll 

0.70000 

1.42857 

2 

MATLAB 

1.05000 

0.31746 

BioCharon 

1.05000 

0.31746 

JigCcll 

1.05000 

0.31746 

3 

MATLAB 

1.16667 

0.17143 

BioCharon 

1.16667 

0.17143 

JigCcll 

1.16667 

0.17143 

4-4  Lysis-Lysogeny  Pathway 

Depending  on  which  operator  site  the  repressor  (A")  binds  to,  the  mutant  phage 
(discussed  in  section  3.7.1)  is  either  virulent  or  temperate.  When  X  binds  to  OR2, 
represented  by  DX-2,  the  synthesis  of  repressor  is  enhanced.  If  X  binds  to  OR3, 
represented  by  DX2,  the  synthesis  of  X  is  inhibited.  High  concentration  levels  of  X 
induce  a  temperate  phage.  On  the  other  hand,  low  concentration  levels  of  repressor 
induce  a  virulent  phage  [10:271]. 

4-4-1  General  Results.  From  equation  3.61,  it  is  obvious  that  the  lysis- 
lysogeny  pathway  can  have  multiple  critical  points.  Three  real  and  two  imaginary 
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critical  points  are  calculated  by  substituting  the  appropriate  values  listed  in  table 
4.3  into  equation  3.61.  Table  4.3  defines  the  following  parameter  values  for  the 
lysis-lysogeny  pathway  system:  kinetic  rate  constants  (k\,  kt,  r,  etc.),  conservation 
constant  (dr),  constant  for  polymerization  of  species  X  (n) ,  and  constant  concen¬ 
tration  for  RNA  polymerase  (P).  Three  real  critical  points,  using  equations  3.58, 
3.59,  and  3.61,  arc  listed  in  table  4.4. 


Table  4.3.  Model  Parameters  (Lysis-Lysogeny  Pathway) 


Parameter 

Value 

k\ 

1 

k- 1 

1 

k2 

1 

k—2 

1 

h 

1 

h 

1 

k_  4 

1 

kt 

0.8 

kd 

6 

r 

0.4 

dj' 

5 

n 

4 

P  =  Po 

1.25 

Table  4.4.  Critical  Points  (Lysis-Lysogeny  Pathway) 


Species 

Critical  Point  1 

Critical  Point  2 

Critical  Point  3 

V 

0.09820 

0.24592 

0.88845 

X‘2 

0.00964 

0.06048 

0.78934 

D 

4.90495 

4.44597 

1.56165 

dx2 

0.04730 

0.26889 

1.23268 

DX* 

0.04730 

0.26889 

1.23268 

dx2x2 

0.00046 

0.01626 

0.97300 

Using  the  same  approach  as  previously  discussed  in  section  3.6,  I  perform  sta¬ 
bility  analysis  on  the  full  and  reduced  systems  of  the  lysis-lysogeny  pathway  (about 
each  critical  point  listed  in  table  4.4).  The  full  system  is  represented  by  equations 
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Table  4.5.  Eigenvalues  for  Full  Lysis-Lysogeny  Pathway 


Critical  Point  1 

Critical  Point  2 

Critical  Point  3 

Eigenvaluei 

0 

0 

0 

Eigenvalue2 

-0.02752 

0.01455 

-0.11060 

Eigenvalue 3 

-0.89098 

-0.83369 

-0.91172 

Eigenvalue 4 

-1.20181 

-1.70225 

-5.02094 

Eigenvalue 5 

-6.43729 

-6.99262 

-7.88717 

Eigenvalue 6 

-11.94908 

-12.12940 

-14.43544 

Table  4.6.  Eigenvalues  for  Reduced  Lysis- Lysogeny  Pathway 


Critical  Point  1 

Critical  Point  2 

Critical  Point  3 

Eigenvalue  1 

-0.02752 

0.01455 

-0.11060 

Eigenvalue  2 

-0.89098 

-0.83369 

-0.91172 

Eigenvalue 3 

-1.20181 

-1.70225 

-5.02094 

Eigenvalue 4 

-6.43729 

-6.99262 

-7.88717 

Eigenvalue 5 

-11.94908 

-12.12940 

-14.43544 

3.50  through  3.55.  The  reduced  system  is  represented  by  equations  3.50,  3.56,  3.52, 
3.57,  and  3.54.  The  calculated  eigenvalues  from  each  stability  analysis  are  listed  in 
tables  4.5  and  4.6. 

With  the  exception  of  the  zero  entries  (Eigenvalue i)  for  each  critical  point 
listed  in  table  4.5,  the  eigenvalues  for  tables  4.5  and  4.6  are  identical.  The  zero 
eigenvalue  listed  for  each  critical  point  in  tabic  4.5  is  expected  because  there  exists  a 
conservation  relation  between  equations  3.52  through  3.55  (i.e. ,  the  full  system  can 
be  reduced).  From  table  4.6,  the  stability  about  (near)  each  critical  point  can  be 
ascertained  (recall:  in  the  reduced  system,  [DX2X2]  =  dr  —  ([D\  +  [DX2\  +  [DX%])). 

As  a  result,  critical  points  one  and  three  are  asymptotically  stable.  This  is 
true  because  all  calculated  eigenvalues  about  critical  points  one  and  three  are  neg¬ 
ative  [7:151].  In  contrast,  critical  point  two  is  not  asymptotically  stable,  because 
Eigenvaluei  is  positive.  As  a  reminder,  small  perturbations  away  from  either  stable 
critical  point  result  in  a  trajectory  back  to  the  respective  critical  point. 
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Conservation  Relation  for  Lysis-Lysogeny  Pathway 


Figure  4.5.  This  illustrates  the  conservation  relation:  dx  =  ([D]  +  [DX 2]  +  [DX £]  + 
[DX2X2]).  It  is  apparent  that  the  sum  of  concentration  levels  among 
these  species  remains  constant  as  the  system  evolves.  In  this  case, 
dx  =  5.  Parameter  settings  and  initial  conditions  are  used  from  tables 
4.3  and  4.7  (trial  4). 

However,  small  perturbations  away  from  the  unstable  critical  point  result  in  a  tra¬ 
jectory  to  one  of  the  stable  critical  points. 

The  existence  of  a  conservation  relation  among  equations  3.52  through  3.55 
implies  that  even  though  individual  concentrations  of  D,  DX 2,  DX£,  and  DX 2X2 
may  change  over  time,  the  sum  of  their  concentrations  remain  constant.  Figure  4.5 
illustrates  this  conservation  rule.  The  plot  (parameter  values  and  initial  conditions 
are  from  tables  4.3  and  4.7  (trial  4),  illustrated  in  figure  4.5,  exemplifies  how  the 
total  sum  of  concentrations  among  species  D.  DX2,  DX2.  and  DX>X>  remain  at  a 
constant  value  (dx)  as  the  modified  Brusselator  system  evolves. 
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4- 4- %  Simulation  Results.  I  perform  a  set  of  four  simulations  (trials)  using 
MATLAB,  BioCharon,  and  JigCell.  The  parameter  values  listed  in  table  4.3  are 
used  in  each  trial.  Initial  conditions  for  each  trial  are  listed  in  table  4.7.  Initial 
species  concentrations  in  trial  1  represent  a  small  perturbation  away  from  critical 
point  1.  Initial  species  concentrations  in  trials  2  and  3  represent  small  perturbations 
away  from  critical  point  2.  Initial  species  concentrations  in  trial  4  represent  a  small 
perturbation  away  from  critical  point  3. 

Simulation  results  from  MATLAB,  BioCharon,  and  JigCell  are  listed  in  table 
4.8.  These  results  are  critical  point  approximations.  Critical  point  approximations 
for  trials  1  and  4  are  similar  to  critical  points  1  and  3  (listed  in  table  4.4).  This 
is  expected  because  trials  1  and  4  are  small  perturbations  away  from  stable  critical 
points. 

Critical  point  approximations  for  trials  2  and  4  are  similar  to  critical  points  1 
and  3,  respectively  (listed  in  table  4.4).  It  is  not  surprising  that  both  trajectories 
move  away  from  critical  point  2,  because  both  trials  represent  small  perturbations 
away  from  an  unstable  critical  point.  In  trial  2,  the  trajectory  path  is  to  critical 
point  1.  In  trial  3,  the  trajectory  path  is  to  critical  point  3.  Apparently,  initial 
conditions  in  trial  2  and  3  result  in  displacements  in  the  attraction  fields  of  critical 
point  1  and  3,  respectively. 

Absolute  error  in  table  4.9  is  calculated  from  entries  listed  in  table  4.8.  Sim¬ 
ulation  results  in  table  4.8  are  critical  point  approximations  from  all  three  software 
tools.  For  both  tables,  M  =  MATLAB,  B  =  BioCharon,  J  =  JigCell.  Absolute  error 
( Ei — defined  in  section  4.2)  for  each  simulation  result  is  listed  in  table  4.9  (columns 
3-8).  Sum  of  absolute  error  {Total e  =  E\  +  E2  +  E%  +  E±  +  E5  +  E§)  is  listed  in 
column  9. 

With  the  exception  of  trial  1,  Total e  for  BioCharon  is  greater  than  Total e  for 
JigCell  in  all  trials.  Maximum  Total e  for  BioCharon  and  JigCell  occur  during  trial 
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2.  The  difference  of  Total e  between  BioCharon  and  JigCcll  is  0.03750;  however, 
Total e  for  BioCharon  is  not  significant  from  a  biological  standpoint. 


Table  4.7.  Initial  Conditions  (Lysis-Lysogeny  Pathway) 


Trial 

X 

AC 

D 

dx2 

DX* 

DX2X2 

1 

0.2000 

0.1000 

4.800 

0.0300 

0.0400 

0.1300 

2 

0.3000 

0.0700 

4.2400 

0.3700 

0.3500 

0.0400 

3 

0.1400 

0.0100 

4.4459 

0.2688 

0.2688 

0.0165 

4 

0.9800 

0.6800 

1.4000 

1.0000 

1.1000 

1.5000 

Table  4.8.  Comparison  Results  (Lysis-Lysogeny  Pathway) 


Trial 

App 

X 

AC 

D 

dx2 

DX* 

dx2x2 

1 

M 

0.09820 

0.00964 

4.90495 

0.04730 

0.04730 

0.00046 

B 

0.09820 

0.00964 

4.90495 

0.04730 

0.04730 

0.00046 

J 

0.09820 

0.00964 

4.90495 

0.04730 

0.04730 

0.00046 

2 

M 

0.88845 

0.78934 

1.56165 

1.23268 

1.23268 

0.97300 

B 

0.88496 

0.78316 

1.56730 

1.22744 

1.22744 

0.96128 

J 

0.88845 

0.78934 

1.56165 

1.23267 

1.23267 

0.97300 

3 

M 

0.09820 

0.00964 

4.90494 

0.04730 

0.04730 

0.00046 

B 

0.09820 

0.00964 

4.90515 

0.04731 

0.04731 

0.00046 

J 

0.09820 

0.00964 

4.90494 

0.04730 

0.04730 

0.00046 

4 

M 

0.88845 

0.78934 

1.56164 

1.23268 

1.23268 

0.97300 

B 

0.88940 

0.79104 

1.56011 

1.23410 

1.23410 

0.97622 

J 

0.88845 

0.78934 

1.56165 

1.23268 

1.23268 

0.97300 

4-5  V.  ftscheri 

4-5.1  General  Results.  Luminescence  in  V.  fischeri  (represented  by  the 
concentration  of  LuxA/B)  is  controlled  by  its  quorum  sensing  system.  This  sys¬ 
tem  has  the  capability  to  ascertain  local  population  density  of  like  bacteria.  When 
the  population  reaches  a  certain  density  level  (i.e.,  a  quorum  of  bacteria),  the  lux 
genes — responsible  for  luminescence — are  subsequently  activated.  In  this  manner, 
luminescence  in  V.  fischeri  exhibits  switch-like  behavior.  Figure  4.6  illustrates  how 
luminescence  in  V.  fischeri  is  activated  when  a  quorum  is  reached. 
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Table  4.9.  Absolute  Error  (Lysis-Lysogeny  Pathway) 


Trial 

App 

E\ 

E2 

e3 

£4 

E5 

Eq 

T  otcilp; 

2 

B 

0.00349 

0.00618 

0.00565 

0.00524 

0.00524 

0.01172 

0.03752 

J 

0.00000 

0.00000 

0.00000 

0.00001 

0.00001 

0.00000 

0.00002 

3 

B 

0.00000 

0.00000 

0.00021 

0.00000 

0.00001 

0.00000 

0.00022 

J 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

4 

B 

0.00095 

0.00170 

0.00153 

0.00142 

0.00142 

0.00322 

0.01024 

J 

0.00000 

0.00000 

0.00001 

0.00000 

0.00000 

0.00000 

0.00001 

As  depicted  in  figure  4.6,  luminescence  in  the  cell  is  negligible  between  zero  and 
twelve  hours,  because  population  density  is  low.  A  quorum  is  apparently  reached 
slightly  after  twelve  hours,  because  the  concentration  of  LuxA/B  begins  to  rapidly 
increase.  In  this  way,  the  ‘switch’  that  controls  luminescence  is  activated  when  a 
quorum  is  achieved. 

“Along  with  LuxR  and  LuxI,  cAMP  receptor  protein  (CRP)  plays  an  important 
role  in  controlling  luminescence  [1:8].”  LuxR  binds  to  Ai  that  forms  the  complex  Co. 
This  complex  subsequently  binds  to  the  lux  box  and  provides  positive  feedback  for 
the  synthesis  of  LuxI  and  negative  feedback  for  the  inibition  of  LuxR.  Increasing 
levels  of  LuxI  enhance  the  production  of  Ai.  CRP  also  binds  to  the  lux  box  and 
provides  positive  feedback  for  the  synthesis  of  LuxR  and  negative  feedback  for  the 
inhibition  of  LuxI  (see  figure  3.4). 

Figures  4.7,  4.8,  and  4.9  illustrate  numerical  simulations  on  the  time  interval — 
[0  40]  hours — of  the  quorum  sensing  system  that  controls  luminescence.  Initial  con¬ 
ditions  for  these  simulations  are  listed  in  table  4.10 — for  trials  1,  2,  and  3.  These 
initial  conditions  are  configured  in  such  a  way  that  all  species  concentrations  are 
identical,  with  the  exception  that  the  concentration  for  CRP  (the  factor),  which 
remains  constant,  but  varies  from  trial  to  trial.  For  trials  1-3,  CRP  is  assigned  the 
initial  concentrations  of  0.7,  20,  and  300  nM  (nanomolar),  respectively. 
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Population  Growth 


xIO'3 


Figure  4.6.  The  top  figure  represents  time  evolution  of  scaled  population  for 

V.  fischeri:  V5  =  1.5e  —  15  liter  and  V  =  1.0e-3  liter.  The  bottom 
figure  represents  time  evolution  of  LnxA/B.  The  concentration  of  this 
protein  directly  correlates  to  the  intensity  of  cell  luminescence. 
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The  following  qualitative  characteristics  are  common  for  all  three  simulations, 
with  varying  CRP  concentrations.  At  low  population  densities,  there  is  a  short 
adjustment  period  (i.e.,  rapid  increase)  for  proteins,  complex  Co,  and  Ai.  The  con¬ 
centrations  for  these  species  remain  nearly  constant  from  3  to  12  hours.  Since  the 
population  density  is  low,  luminescence  is  low  (indicated  by  LuxA/B).  Shortly  after 
12  hours  (for  trials  2  and  3),  another  adjustment  is  initiated,  where  there  is  another 
rapid  increase  in  proteins,  Co,  and  autoinducer  Ai.  After  20  hours,  the  concentration 
for  these  species  approach  nearly  constant  levels.  Luminescence  in  the  cell  is  at  its 
maximum  level. 

Figures  4.7,  4.8,  and  4.9  also  illustrate  how  varying  concentrations  of  CRP  (i.e., 
CRP  at  low,  moderate,  and  high  concentrations)  possibly  affect  luminescence  in  V. 
hscheri.  CRP  binds  to  the  lux  box  and  provides  positive  feedback  to  Ol  (luxR)  and 
negative  feedback  to  Or  (luxICDABEG).  The  concentration  of  CRP  determines  the 
degree  of  positive  and  negative  feedback  asserted  to  Ol  and  Or. 

When  the  concentration  of  CRP  equals  20,  increasing  Co  provides  positive 
feedback  to  Or.  This  causes  the  production  rate  of  LuxI  and  Ai  to  increase.  As  the 
concentration  levels  of  these  species  increase,  the  total  luminescence  increases  (i.e., 
LuxA/B  increases).  In  this  case,  cell  luminescence  is  high.  Figure  4.8  illustrates  this 
behavior. 

When  the  concentration  of  CRP  equals  300,  a  high  degree  of  negative  feed¬ 
back  is  applied  to  Or.  This  negative  feedback  delays  the  synthesis  of  LuxI  and  Ai. 
Therefore,  overall  luminescence  is  moderate.  Figure  4.9  illustrates  this  behavior. 

When  the  concentration  of  CRP  equals  0.7,  positive  feedback  to  Or  is  low.  As 
a  result,  a  significantly  lower  concentration  of  LuxR  is  initially  produced.  Conse¬ 
quently,  the  concentration  level  of  Co  is  not  sufficient  to  activate  Or.  In  this  case, 
cell  luminescence  is  extremely  low.  Figure  4.7  illustrates  this  behavior. 
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Concentration  [nM] 


CRP  =  0.7  nM 


[hours] 

Figure  4.7.  This  illustrates  a  numerical  simulation  (log- linear  plot)  for  LuxR, 
LnxA/B,  LuxI,  Ai,  and  Co  in  the  bioluminescence  control  model  for 
V.  fischeri.  In  this  case,  CRP  =  0.7  (trial  1). 
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Concentration  [nM] 


4-5.2  Comparison  Results.  I  compare  simulation  results — for  the  biolu¬ 
minescence  control  model — generated  by  MATLAB  and  JigCcll.  The  metrics  are 
absolute  error  and  relative  error.  All  species  concentrations  are  in  nano  molar  con¬ 
centrations  (nM),  except  for  population.  The  unit  of  measure  for  population  is  the 
scaled  volume  of  bacteria  (populationxvb/V ,  where  ty  =  1.5E-15  liter  and  V  =  1.0E-3 
liter  [4:4-5]).  A  set  of  six  simulations  are  performed  on  JigCell  and  MATLAB. 

Initial  concentrations  for  all  species — separated  by  trial — are  listed  in  table 
4.10.  The  initial  concentrations  for  each  trial  are  identical,  with  the  exception  of 
CRP.  For  trials  one  and  four,  the  concentration  of  CRP  is  0.7.  For  trials  two  and 
five,  the  concentration  of  CRP  is  20.  For  trials  three  and  six,  the  concentration  of 
CRP  is  300. 

Simulations  results  performed  on  the  time  interval — [0  5]  hours — are  listed  in 
tables  4.11,  4.12,  and  4.13.  Simulation  results  performed  on  the  time  interval — [20 
25]  hours — are  listed  in  tables  4.14,  4.15,  and  4.16.  With  the  exception  of  trial  5 
(listed  in  table  4.15),  the  largest  absolute  error  between  MATLAB  and  JigCell  occurs 
for  species  LuxR.  During  trial  5,  the  largest  absolute  error  occurs  for  species  Ai.  This 
is  not  surprising  because  Ai  is  significantly  larger  than  all  other  species. 

On  the  other  hand,  the  relative  error  is  minimal  for  all  trials.  The  maximum 
relative  error  occurs  during  trials  1  and  3  for  scaled  population  density.  Nonetheless, 
the  maximum  error  is  only  0.559%. 

4-6  Glycolytic  Pathway 

4-6.1  General  Results.  As  discussed  in  section  3.7.4,  each  of  the  10  reac¬ 
tions  or  steps  in  the  glycolytic  pathway  (illustrated  in  figure  A)  is  catalyzed  by  a 
specific  enzyme.  In  order  for  a  reaction  to  take  place  at  each  step,  a  specified  enzyme 
must  be  present  to  interact  with  the  substrate.  If  the  required  enzyme  is  absent  at  a 
certain  step,  that  reaction  is  not  able  to  take  place.  I  incorporate  this  idea  into  the 
design  of  experimental  trials  conducted  in  this  section. 
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Table  4.10.  Initial  Conditions  (V.  fischeri) 


Species 

Trial  1 

Trial  2 

Trial  3 

Trial  4 

Trial  5 

Trial  6 

population 

1.5E-07 

1.5E-07 

1.5E-07 

0.00149 

0.00149 

0.00149 

luxR 

0 

0 

0 

0.98132 

0.37028 

2.58986 

luxICDABEG 

0 

0 

0 

0.00063 

4.81455 

0.39765 

LuxR 

0 

0 

0 

869.16824 

265.33072 

4241.58358 

LuxI 

0 

0 

0 

0.54098 

4242.81014 

285.69095 

LuxA/B 

0 

0 

0 

0.54098 

4242.81014 

285.69095 

LuxC/D/E 

0 

0 

0 

0.09353 

720.96687 

58.80923 

Ai 

6 

6 

6 

1.38029 

3957.78272 

104.74681 

CO 

0 

0 

0 

0.11666 

102.06949 

43.03559 

CRP 

0.7 

20 

300 

0.7 

20 

300 

Table  4.11.  Trial  1  Simulation  Results  (V.  fischeri) 


Species 

MATLAB 

JigCell 

absolute  error 

relative  error 

population 

4.80855E-06 

4.78166E-06 

0.00000 

0.00559 

luxR 

0.97021 

0.97021 

0.00000 

0.00000 

luxICDABEG 

0.00127 

0.00127 

0.00000 

0.00000 

LuxR 

515.95140 

515.93890 

0.01250 

0.00002 

LuxI 

0.75306 

0.75297 

0.00009 

0.00012 

LuxA/B 

0.75306 

0.75297 

0.00009 

0.00012 

LuxC/D/E 

0.17575 

0.17577 

0.00002 

0.00011 

Ai 

3.64003 

3.64009 

0.00006 

0.00002 

CO 

0.17985 

0.17985 

0.00000 

0.00000 

Table  4.12.  Trial  2  Simulation  Results  (V.  fischeri) 


Species 

MATLAB 

JigCell 

absolute  error 

relative  error 

population 

4.80875E-06 

4.78166E-06 

0.00000 

0.00563 

luxR 

9.80402 

9.80402 

0.00000 

0.00000 

luxICDABEG 

0.04130 

0.04130 

0.00000 

0.00000 

LuxR 

5208.60414 

5208.48680 

0.11734 

0.00002 

LuxI 

24.55205 

24.54652 

0.00553 

0.00023 

LuxA/B 

24.55205 

24.54652 

0.00553 

0.00023 

LuxC/D/E 

5.71036 

5.7101 

0.00022 

0.00004 

Ai 

3.66911 

3.66896 

0.00015 

0.00004 

CO 

1.82981 

1.82969 

0.00012 

0.00007 
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Table  4.13.  Trial  3  Simulation  Results  (V.  fischeri) 


Species 

MATLAB 

JigCell 

absolute  error 

relative  error 

population 

4.81013E-06 

4.78166E-06 

0.00000 

0.00592 

luxR 

14.10938 

14.10937 

0.00001 

0.00000 

luxICDABEG 

0.00820 

0.00820 

0.00000 

0.00000 

LuxR 

7488.10943 

7487.93700 

0.17243 

0.00002 

LuxI 

4.87340 

4.87261 

0.00079 

0.00016 

LuxA/B 

4.87340 

4.87261 

0.00079 

0.00016 

LuxC/D/E 

1.13434 

1.13441 

0.00007 

0.00006 

Ai 

3.64504 

3.64506 

0.00002 

0.00001 

CO 

2.61356 

2.61350 

0.00006 

0.00002 

Table  4.14.  Trial  4  Simulation  Results  (V.  fischeri) 


Species 

MATLAB 

JigCell 

absolute  error 

relative  error 

population 

0.00150 

0.00150 

0.00000 

0.00000 

luxR 

0.98142 

0.98142 

0.00000 

0.00000 

luxICDABEG 

0.00070 

0.00070 

0.00000 

0.00000 

LuxR 

882.73431 

882.72736 

0.00695 

0.00001 

LuxI 

0.61346 

0.61345 

0.00001 

0.00002 

LuxA/B 

0.61346 

0.61345 

0.00001 

0.00002 

LuxC/D/E 

0.10408 

0.10408 

0.00000 

0.00000 

Ai 

1.45469 

1.45469 

0.00000 

0.00000 

CO 

0.12489 

0.12489 

0.00000 

0.00000 

Table  4.15.  Trial  5  Simulation  Results  (V.  fischeri) 


Species 

MATLAB 

JigCell 

absolute  error 

relative  error 

population 

0.00150 

0.00150 

0.00000 

0.00000 

luxR 

0.30227 

0.30224 

0.00003 

0.00010 

luxICDABEG 

4.84908 

4.84909 

0.00001 

0.00000 

LuxR 

166.69868 

166.70787 

0.00919 

0.00006 

LuxI 

4358.01106 

4357.98780 

0.02326 

0.00001 

LuxA/B 

4358.01106 

4357.98780 

0.02326 

0.00001 

LuxC/D/E 

727.25376 

727.25336 

0.00040 

0.00000 

Ai 

6989.94902 

6989.86470 

0.08432 

0.00001 

CO 

113.33358 

113.34041 

0.00683 

0.00006 
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Table  4.16.  Trial  6  Simulation  Results  (V.  fischeri) 


Species 

MATLAB 

JigCell 

absolute  error 

relative  error 

population 

0.00150 

0.00150 

0.00000 

0.00000 

luxR 

1.45400 

1.45381 

0.00019 

0.00013 

luxICDABEG 

0.43556 

0.43556 

0.00000 

0.00000 

LuxR 

1378.32756 

1378.45730 

0.12974 

0.00009 

LuxI 

387.46805 

387.46323 

0.00482 

0.00001 

LuxA/B 

387.46805 

387.46323 

0.00482 

0.00001 

LuxC/D/E 

65.24738 

65.24728 

0.00010 

0.00000 

Ai 

447.74904 

447.74402 

0.00502 

0.00001 

CO 

59.98567 

59.99049 

0.00482 

0.00008 

In  order  to  present  nontrivial  results,  pyruvate  kinase  (1/35)  is  ‘removed’  in 
order  to  nullify  reaction  ten.  If  reaction  ten — a  one-way  reaction — is  permitted  to 
take  place,  the  concentrations  of  all  species  located  upstream  from  reaction  ten  would 
be  eventually  depleted.  Since  intermediate  reactions  five,  six,  seven,  and  nine  are 
two-way  reactions,  it  is  expected  that  those  species  will  approach  non- zero  values  as 
the  system  evolves  towards  a  critical  point.  Therefore,  I  implement  a  computational 
model  that  contains  the  first  34  equations  listed  in  equation  3.68. 

As  discussed  in  section  4.1,  the  computational  cost  to  obtain  accurate  simula¬ 
tion  results  from  MATLAB  is  substantially  less  than  from  BioCharon  and  JigCcll. 
Therefore,  simulations  for  glycolytic  pathway  are  performed  only  using  MATLAB. 

BioCharon  generates  the  same  system  of  ODEs  (computational  model)  listed 
as  equation  3.68.  However,  JigCcll  generates  a  reduced  computational  model,  con¬ 
sisting  of  16  dependent  species  and  18  independent  species  (state  variables).  This 
occurs  because  JigCell  automatically  identifies  conservation  relations  among  species. 
Conservation  relations  automatically  generated  by  JigCell  are  listed  in  equation  4.1 
(T,  is  a  constant  determined  by  the  initial  conditions  of  corresponding  species). 

Figure  4.10  illustrates  the  conservation  relation  listed  in  equation  4.1  (expres¬ 
sion  in  which  r/13  is  the  dependent  variable).  The  concentrations  of  the  relevant 
species  at  each  time  step  are  obtained  by  performing  a  simulation  of  the  glycolytic 
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Conservation  Relation  in  Glycolytic  Pathway 
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Figure  4.10.  This  illustrates  the  conservation  relation  identified  by  JigCell,  listed 
in  equation  4.1.  This  is  the  expression  in  which  y18  is  the  dependent 
variable.  In  this  case  the  algebraic  expression  of  evolving  species 
remains  constant  at  -2000  (identified  as  T13  in  table  4.18). 


pathway,  using  initial  concentrations  listed  in  table  4.17.  This  plot  validates  that 
JigCell  correctly  identified  a  relatively  complicated  conservation  relation,  consist¬ 
ing  of  16  species.  In  this  case,  the  conservation  constant  among  those  species  is 
T13  =  —2000.  All  relation  constants  are  listed  in  table  4.18.  The  values  for  these 
constants  are  derived  from  the  initial  species  concentrations  substituted  into  the 
appropriate  expression  in  equation  4.1. 


Vi  =  (T9  +  0.5  *y2  +  0.5  *y5  +  0.5  *y8~  0.5  *  yw  -  0.5  *  y28- 
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0.5  *  1/30  +  0.5  *  1/32  -  0.5  *  1/34  -  0.5  *y7  -  0.5  *  y31) 

2/6  =  (^16  —  2/5) 

2/9  =  (T10  -  2/s) 

2/12  =  (^4  —  2/ll) 

2/13  =  (T13  -  0.5  *  2/2  +  0.5  *  2/5  +  0.5  *  y8  -  0.5  *  y10  -  y12  +  yu~ 

0.5  *  yw  +  0.5  *  2/is  +  0.5  *  y20  +  0.5  *  y2 s  +  0.5  *  //16+ 

0-5  *  2/16  +  0.5  *  1/34  —  0.5  *  2/7  +  0.5  *  y31  —  0.5  *  y17) 

2/15  —  (^5  —  2/m) 

2/19  =  (^6  —  //is) 

2/21  =  (Ti5  +  2/20) 

2/22  =  (7\  —  1/20  +  2/23) 

2/24  =  (^7  —  2/23) 

2/25  =  (Til  +  2/2  -  2/5  +  2/3  +  2/12  -  2/20  +  2/23  — 

2/28  —  2/30  +  2/32  —  2/34  —  2/3l) 

2/26  =  (^2  +  2/2  —  2/5  +  2/3  +  2/12) 

2/27  —  (^8  —  2/2  +  2/5  —  2/3  —  2/12) 

2/30  =  (^14  —  2/29) 

2/32  =  “(T3  +  2/2  +  2/3  -  2/4) 

2/33  =  (^12  —  2/32)  (4.1) 

4-6.2  Comparison  Results.  I  compare  simulation  results  between  the  com¬ 
putational  model  generated  by  JigCcll  (consisting  of  18  state  variables)  and  the 
computational  model  listed  as  equation  3.68  (consisting  of  34  state  variables).  Ta¬ 
ble  4.17  lists  initial  species  concentrations  for  simulation  results  listed  in  table  4.19. 
Simulation  results  from  the  full  and  reduced  models  are  listed  in  columns  one  and 
two  of  table  4.19.  Absolute  error  is  listed  in  column  three.  The  unit  of  measure  for 
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species  concentrations  can  be  in  any  consistent  unit.  The  simulation  time  is  1000 
seconds. 

The  maximum  absolute  error  occurs  for  species  t/2  (listed  in  table  4.19).  This 
error  is  negligible,  because  the  error  is  just  0.0003%.  Therefore,  the  reduced  model 
is  validated  by  MATLAB  results. 


Table  4.17.  Initial  Conditions  (Glycolytic  Pathway) 


Species 

[Value] 

Species 

[Value] 

Species 

[Value] 

Species 

[Value] 

U\ 

500 

2/2 

1000 

2/3 

1000 

2/4 

1000 

Vs 

1000 

2/6 

0.00001 

2/7 

0.00001 

2/8 

1000 

2/9 

0.00001 

2/io 

0.00001 

2/ii 

1000 

2/12 

0.00001 

2/13 

0.00001 

2/14 

1000 

2/15 

0.00001 

2/16 

0.00001 

2/17 

0.00001 

2/18 

1000 

2/19 

0.00001 

2/20 

1000 

2/21 

1000 

2/22 

0.00001 

2/23 

1000 

2/24 

0.00001 

2/25 

0.00001 

2/26 

1000 

2/27 

0.00001 

2/28 

0.00001 

2/29 

1000 

2/30 

0.00001 

2/31 

0.00001 

2/32 

1000 

2/33 

0.00001 

2/34 

0.00001 

4-7  Yeast  Cell  Cycle 

4-7.1  General  Results.  As  discussed  in  section  3.7.3,  the  cell  division 
cycle  of  budding  yeast  is  characterized  by  the  system  alternating  between  two  self- 
maintaining  states — State  1  and  State2-  CDK  activity  (combined  activities  of  Clb’s 
and  Cln’s)  is  low  in  State  1,  while  their  enemies — Sicl  and  Hctl  are  high.  Therefore, 
concentration  levels  of  the  Clb’s  and  Cln’s  remain  mostly  low  during  State  1.  On 
the  other  hand,  CDK  activity  is  high  in  State 2.  As  a  result,  the  concentrations  of 
the  Clb’s  and  Cln’s  are  high,  where  their  enemies’  concentrations  remain  mostly  low 
during  State 2. 

Control  is  transferred  from  State  1  to  State 2  (and  vice  versa)  by  means  of  two 
transitions — Start  and  Finish.  In  the  molecular  control  model  of  CDK  activities 
(figure  3.5),  Start  is  initiated  when  the  concentration  of  species  ORI  reaches  a  con¬ 
centration  of  1  unit  (approached  from  below)  [1:374],  Finish  is  initiated  when  species 
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Tabic  4.18. 


Conservation  Relation  Constants  (Glycolytic  Pathway) 


Constant 

Value 

T\ 

0.00001 

t2 

-0.00001 

t3 

-2000 

t4 

1000.000010 

n 

1000.000010 

n 

1000.000010 

T7 

1000.00001 

t8 

1000.000020 

t9 

-1499.999970 

Tw 

1000.000010 

T\\ 

-1999.999960 

Ti2 

1000.000010 

Tl3 

-1999.999980 

Tu 

1000.000010 

Ti5 

0 

Tin 

1000.000010 

SPN  reaches  a  concentration  of  1  unit  (approached  from  below)  [1:374],  Simulation 
results,  generated  by  MATLAB  from  equations  3.66  and  3.67,  depict  important  char¬ 
acteristics  of  molecular  interaction  in  the  budding  yeast  cell  division  cycle.  These 
results  are  illustrated  in  figures  4.11,  4.12,  and  4.13. 

Activities  of  Clb2T,  Clb5T,  and  Cln2  are  illustrated  in  figure  4.11.  Collectively, 
the  activities  of  these  species  are  low  during  State i  and  high  during  State 2.  Activities 
of  SiclT  and  HiclT  (enemies  of  CDKs)  are  illustrated  in  figure  4.12.  Their  activities 
are  opposite  from  the  CDKs.  Figure  4.13  illustrates  cell  division.  Cell  division  occurs 
around  380  minutes.  Rules  for  cell  division  are  listed  in  the  caption  of  figure  4.13. 

4-7.2  Comparison  Results.  In  the  previous  sections,  initial  species  concen¬ 
trations  are  perturbed  during  each  trial.  In  this  section,  the  reaction  rate  constant, 
kd, 20,  is  perturbed  from  a  value  of  0.080  to  0.092  at  0.006  increments.  This  rate 
constant  is  listed  in  equation  3.66  (in  the  differential  equation  defining  the  rate  of 
change  for  species  Cdc20).  Initial  species  concentrations  and  other  model  pararne- 
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Table  4.19.  Comparison  Results  (Glycolytic  Pathway) 


Species 

Full  Model 

Reduced  Model 

Absolute  Error 

1.0E+003  * 

1.0E+003  * 

1.0E-005  * 

V\ 

0.00000000000000 

0.00000000000001 

0.00000112688929 

V2 

0.32719541332651 

0.32719541606607 

0.27395537358643 

V3 

1.01346976779097 

1.01346976794359 

0.01526191226731 

Vi 

0.01673690254825 

0.01673690290803 

0.03597766031760 

Vs 

1.00000001000000 

1.00000000999480 

0.00052004907047 

Ve 

0.00000000000000 

0.00000000000520 

0.00052003770179 

Vi 

0.00000000000000 

0.00000000000001 

0.00000067220655 

Vs 

1.00000001000000 

1.00000000999543 

0.00045685055738 

V9 

0.00000000000000 

0.00000000000457 

0.00045682782002 

yio 

0.00000000000000 

0.00000000000000 

0.00000000010772 

2/n 

1.00000001000000 

1.00000000999601 

0.00039946144170 

2/ 12 

0.00000000000000 

0.00000000000399 

0.00039934775486 

2/13 

0.00000000013785 

0.00000000013526 

0.00025951412943 

2/14 

0.99999959644638 

0.99999959644794 

0.00015545538190 

2/15 

0.00000041355362 

0.00000041355206 

0.00015550629033 

2/16 

0.00000334154936 

0.00000334154891 

0.00004481718911 

2/17 

0.00000037128326 

0.00000037128321 

0.00000510898911 

2/18 

0.99888739948217 

0.99888739963078 

0.01486117753302 

2/19 

0.00111261051783 

0.00111261036922 

0.01486108365256 

2/20 

0.00111698053804 

0.00111698060354 

0.00654993137505 

2/21 

0.00111698053804 

0.00111698060354 

0.00654993135285 

2/22 

0.98653033025512 

0.98653033007102 

0.01840957111199 

2/23 

0.98764730047205 

0.98764730067456 

0.02025145704465 

2/24 

0.01235270952795 

0.01235270932544 

0.02025145473539 

2/25 

0.00000003803393 

0.00000003803361 

0.00003179414129 

2/26 

0.34066516520734 

0.34066516401886 

0.11884819173247 

2/27 

0.65933484479266 

0.65933484598114 

0.11884818604813 

2/28 

0.00000106027192 

0.00000106027190 

0.00000260618866 

2/29 

0.99682927970533 

0.99682927978787 

0.00825422148409 

2/30 

0.00317073029467 

0.00317073021213 

0.00825437611596 

2/31 

0.00000954244732 

0.00000954244706 

0.00002542509883 

2/32 

0.67607171776758 

0.67607171889837 

0.11307893146295 

2/33 

0.32392829223242 

0.32392829110163 

0.11307886893519 

2/34 

0.00008588202584 

0.00008588202355 

0.00022895867979 

concentration 


Q  I _ Lfcl _ 1 _ ~~l  — I - - C  1  - 1 - Li--- _ I 

280  300  320  340  300  380  400  420  440  460  480 


minutes 

Figure  4.11.  This  depicts  the  activitiy  levels  of  Gib’s  and  Cln’s.  During  State i, 
the  activity  of  both  cyclins  are  low.  However,  their  activity  is  high 
during  State 2. 
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280  300  320  340  360  380  400  420  440  460 

minutes 

Figure  4.12.  This  depicts  the  activity  levels  of  Sicl  and  HCT1 — enemies  of  Gib’s. 

Their  activity  levels  are  opposite  of  the  Gib’s  and  Gin’s,  displayed  in 
figure  4.11. 
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concentration 
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minutes 

Figure  4.13.  Cell  division  occurs  around  380  minutes.  Cell  division  is  initiated 
when  the  concentration  of  Clb2  drops  below  a  threshold  of  0.3  con¬ 
centration  units.  At  cell  division,  the  mass  is  divided  between  the 
mother  and  daughter  cell  as  follows: 

(mass  of  daughter  cell  at  birth)  =  0.433  x  (mass  at  cell  separation) 
and  (mass  of  mother  cell  at  birth)  =  (1  -  0.433)  x  (mass  at  cell 
separation)  [8:374], 
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ters  remain  constant  for  all  three  trials  (values  are  listed  in  appendix  A — xpp  output 
file). 

The  metrics  for  comparison  results  are  absolute  error  and  sum  of  absolute  error 
(between  MATLAB  and  JigCell  simulation  results).  “All  concentration  variables  are 
scaled  so  that  their  maximal  values  are  pure  numbers  of  order  1  [1:386].”  The 

sampled  time  interval  is  [280  480]  minutes. 

Tables  4.20  and  4.21  list  the  percentage  of  change  caused  by  perturbing  k^2 o- 
In  both  tables  listing  simulation  results  generated  by  JigCell  and  MATLAB,  respec¬ 
tively,  SiclT  is  the  most  sensitive  parameter.  SiclT  is  decreased  -18.1%  on  JigCell 
and  decreased  -17.1%  in  MATLAB. 

Comparison  results  (between  JigCell  and  MATLAB)  are  listed  in  table  4.22. 
The  length  parameter  is  the  interval  immediately  after  cell  division  and  right  before 
cell  division  occurs.  For  all  trials,  maximum  absolute  error  occurs  in  Length.  The 
maximum  sum  of  absolute  error  occurs  during  trial  3. 


Table  4.20.  MATLAB  kd.20  Perturbation  Results  (Budding  Yeast) 


Name 

Trial  1 

Trial  2 

Trial  3 

Overall  Change 

cmT 

1.48672 

1.52294 

1.54121 

3.7% 

ClbhT 

0.26888 

0.26682 

0.26644 

-0.9% 

Cln2 

0.52380 

0.51813 

0.51442 

-1.8% 

Siclx 

1.40809 

1.29039 

1.15385 

-18.1% 

Hctl 

0.99995 

0.99995 

0.99989 

0.0% 

Length 

145.50000 

144.10000 

143.79830 

-1.2% 

kd,  20 

0.080 

0.086 

0.092 

15.0% 

4-  8  Summary 

In  this  chapter  I  perform  numerical  simulations  for  the  following  models:  mod¬ 
ified  Brusselator,  lysis-lysogeny  pathway  of  a  mutant  bacteriophage,  quorum  sensing 
system  that  controls  luminescence  in  V.  fischeri,  a  glycolytic  pathway,  and  control 
system  for  CDKs  in  budding  yeast  cell  division  cycle.  I  present  general  results 
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Table  4.21.  JigCell  2o  Perturbation  Results  (Budding  Yeast) 


Name 

Trial  1 

Trial  2 

Trial  3 

Overall  Change 

cmT 

1.48167 

1.50582 

1.52887 

3.2% 

ClbbT 

0.26897 

0.26767 

0.26654 

-0.9% 

Cln2 

0.52345 

0.51816 

0.51367 

-1.9% 

Sict 

1.39631 

1.26793 

1.15715 

-17.1% 

Hctl 

0.99986 

0.99985 

0.99985 

0.0% 

Length 

144.90000 

144.90000 

144.90000 

0.0% 

kd,  20 

0.080 

0.086 

0.092 

15.0% 

Table  4.22.  Comparison  Results  (Budding  Yeast) 


Name 

Trial  1 

Trial  2 

Trial  3 

cmT 

0.00505 

0.01712 

0.01234 

ClbhT 

0.00009 

0.00085 

0.00010 

Cln2 

0.00035 

0.00003 

0.00075 

Sict 

0.01178 

0.02246 

0.00330 

Hctl 

0.00009 

0.00010 

0.00004 

Length 

0.60000 

0.80000 

1.10170 

Total  Error 

0.61735 

0.84055 

1.11824 

that  depict  important  characteristics  of  each  model.  In  order  to  perform  a  quasi- 
performance  test  for  BioCharon  and  Jigcell,  I  compare  simulation  results  generated 
by  these  software  tools. 

The  models  for  V.  hscheri  and  budding  yeast  are  based  on  journal  publications 
(references  [1],  [4],  [8]).  I  successfully  duplicate  the  published  results. 
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V.  Summary  and  Conclusion 


5. 1  Summary 

The  overall  goal  of  this  thesis  was  to  implement  mathematical  models  (via  rate- 
equation  approach)  for  three  well-documented  cellular  systems  using  BioCharon, 
JigCell,  and  MATLAB.  Each  model  characterizes  a  specific  control  mechanism  that 
synchronizes  a  set  of  intracellular  processes  within  a  genetic  regulatory  network.  The 
experience  gained  from  implementing  these  models  is  intended  to  contribute  to  Air 
Force  toxicology  studies  for  the  following  reason.  “Mathematical  models  are  useful 
for  providing  a  framework  for  integrating  data  and  gaining  insights  into  the  static 
and  dynamic  behavior  of  complex  biological  systems  such  as  networks  of  interacting 
genes  [26:247].'"  In  short,  the  augmentation  of  intracellular  models  into  Air  Force 
toxicology  studies  has  the  potential  to  significantly  improve  the  quality  and  efficiency 
of  the  studies. 

In  order  to  achieve  my  overall  goal,  I  completed  the  following  tasks. 

1.  Acquired  a  sound  understanding  of  the  fundamentals  of  biochemistry  and  mi¬ 
crobiology  that  are  related  to  biological  reaction  systems. 

2.  Implemented  computational  (mathematical)  models  via  rate-equation  approach 
for  the  following  molecular  systems: 

•  modified  Brusselator 

•  glycolyic  pathway 

•  lysis- lysogeny  pathway  of  bacteriophage 

•  bioluminescence  control  in  V.  fischeri 

•  cell  cycle  control  in  budding  yeast 

3.  Became  proficient  in  using  three  software  tools:  BioCharon,  JigCell,  and  MAT- 
LAB. 
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In  addition,  simulation  results  generated  by  BioCharon  and  JigCell  were  compared 
to  simulation  results  generated  by  MATLAB  to  gauge  the  accuracy  of  these  two  trial 
releases  by  DARPA. 

5.2  Conclusion 

Reaching  the  overall  goal  of  this  thesis  has  led  to  a  number  of  conclusions.  The 
rate-equation  approach  is  appropriate  for  deterministically  modeling  both  small  and 
large  molecular  systems  (e.g.,  the  modified  Brusselator  and  the  glycolytic  pathway). 
There  is  a  clear  need  for  the  use  of  software  modeling  tools  that  automate  all  or 
some  of  the  steps  in  the  rate-equation  approach.  Without  such  tools,  the  modeler 
easily  becomes  overburdened  in  attempting  to  derive  a  computational  model  from 
other  than  a  simplistic  biological  system. 

From  this  thesis,  the  merits  of  BioCharon  and  JigCell  are  apparent.  The 
graphical  approach — implemented  by  BioCharon — is  highly  intuitive  and  represents 
the  higher  level  of  abstraction  in  building  models.  By  means  of  a  GUI,  the  user  is 
able  to  quickly  create  a  reaction  system  by  simply  manipulating  two  main  objects 
(nodes  and  arcs).  The  spreadsheet  approach — implemented  by  JigCell — is  relatively 
user-friendly;  however,  the  user  must  enter  the  entire  set  of  reaction  equations  to 
characterize  a  reaction  system.  Both  software  tools  allow  the  user  to  view  the  com¬ 
putational  model  generated  from  user  input.  This  is  crucial  in  the  debugging  process. 

The  three  primary  cellular  models — lysis-lysogeny  pathway,  bioluminescence 
control,  and  cell  cycle  control —  are  outstanding  learning  aids  in  comprehending 
gene  expression  (when  a  gene  is  active  during  transcription).  Gene  expression  is 
represented  at  various  levels  of  detail  in  each  model.  For  the  lysis-lysogeny  pathway 
model,  transcription  and  translation  are  combined  into  equation  3.41.  The  left  hand 
side  of  equation  3.41  depicts  transcription  of  mRNA  transcripts,  and  the  right  hand 
side  depicts  translation  of  mRNA  transcripts  into  n  protein  molecules.  For  the 
bioluminescence  control  model,  the  rates  for  synthesis  of  mRNA  transcripts  from 
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the  left  and  right  operons  ( Or  and  Or)  are  depicted  by  x\  and  X2  in  equation 
3.62.  Rates  for  translation  of  messenger  RNA’s  into  proteins  are  depicted  by  x:il 
X4,  X5,  and  Xfj  in  the  same  equation.  For  the  cell  cycle  control  model,  transcription 
and  translation  rates  are  collectively  modulated  by  transcription  factors  (e.g.,  the 
activation  level  of  transcription  factor  SBF  ehances  the  synthesis  of  Cln2  in  figure 
3.5). 

Gene  expression  is  direcly  relevant  to  toxicology  studies.  “Almost  without 
exception,  gene  expression  is  altered  during  toxicity,  as  either  a  direct  or  indirect 
result  of  toxicant  exposure.  The  challenge  facing  toxicologists  is  to  define,  under  a 
given  set  of  experimental  conditions,  the  characteristic  and  specific  pattern  of  gene 
expression  elicited  by  a  given  toxicant  [21:153].”  Since  gene  expression  can  be  char¬ 
acterized  by  mathematical  models,  specific  models  can  be  used  to  give  predictions 
on  how  toxicant  exposures  affect  cells. 

5.3  Recommendations 

Further  research  into  the  construction  and  evaluation  of  intracellular  models 
would  benefit  Air  Force  toxicology  studies.  In  gaining  the  necessary  experience,  five 
reaction  systems  were  thoroughly  studied  and  implemented.  The  next  step  could  be 
to  construct  a  model  that  is  coupled  with  experimental  data  (previously  discussed 
section  2.1).  Such  data  could  be  freely  obtained  from  the  public  domain  and/or  from 
the  Air  Force. 

The  use  of  deterministic  models  to  depict  reaction  systems  is  not  always  ap¬ 
propriate.  “Conventional  deterministic  kinetics  cannot  be  used  to  predict  statistics 
of  regulatory  systems  that  produce  probabilistic  outcomes.  Rather,  a  stochastic 
kinetic  analysis  must  be  used  to  predict  statistics  of  regulatory  outcomes  for  such 
stochastically  regulated  systems  [3:1633].”  Therefore,  it  is  worthwhile  to  thoroughly 
study  and  implement  mathematical  models  using  the  stochastic  kinetics  approach 
(previously  discussed  in  section  2.8). 
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Even  though  BioCharon  and  JigCell  have  user-friendly  front  ends  for  modeling 
reaction  systems,  their  respective  ODE  solvers  need  improvement.  The  computa¬ 
tional  cost  for  both  tools  was  larger  than  that  for  MATLAB  when  attempts  were 
made  to  implement  models  for  the  bioluminescence  control  and  glycolytic  pathway 
models.  Once  an  ODE  solver  that  is  comparable  to  MATLAB’s  accuracy  and  effi¬ 
ciency  is  incorporated  into  both  software  packages,  both  will  become  very  desirable 
tools  for  modeling  and  analyzing  deterministic  systems. 
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Appendix  A.  JigCell  Model  Builder  Entries 

In  this  section,  I  present  my  entries  for  each  model  spreadsheet  (using  JigCell 
Model  Builder)  and  the  corresponding  output  hies  (full  xpp)  that  characterize  the 
following  models:  lysis-lysogeny  pathway,  quorum  sensing  system,  and  a  control 
system  for  the  budding  yeast  cell  cycle.  The  spreadsheets  are  illustrated  in  the 
following  figures. 

Figure  A.l:  lysis-lysogeny  pathway 
Figure  A. 2:  quorum  sensing  system 
Figure  A. 3:  control  system  for  budding  yeast  cell  cycle 

Each  output  hie  is  listed  immediately  after  its  corresponding  model  spreadsheet. 

I  only  list  entries  for  the  Reaction  and  Modifiers  and  Contants  columns  from 
each  spreadsheet.  The  user  enters  the  set  of  reaction  equations  in  the  Reaction 
column.  In  the  Modifiers  and  Contants  column,  the  user  defines  each  rate  law. 

An  important  aspect  of  user  interaction  with  the  model  spreadsheet  is  that 
the  user  cannot  explicitly  define  a  reversible  reaction.  Consequently,  the  user  must 
divide  a  reversible  reaction  into  two  separate,  one-way  reactions.  For  example,  the 
following  two  reaction  equations  are  the  entries  for  D  +  X2  kl2  DX2  (equation  3.38). 

1.  D  +  X2  ->  DX2 

2.  DX2  ->  D  +  X2 

Entries  (user  input  is  underlined)  in  the  Modifiers  and  Contants  column  for  reactions 
1  and  2,  respectively,  are  Kf  =  E  and  Kf  =  k- 2.  The  user  must  assign  numerical 
values  to  k2  and  k-2  in  a  separate  spreadsheet  (the  Constants  spreadsheet).  Detailed 
information  about  JigCell  Model  Builder  can  be  found  in  [30]. 
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Figure  A.l. 


Reaction 

Modifiers  and  Constants 

1 

->X 

Kf=0.4 

2 

2X->X2 

KM 

3 

X2->2X 

Kf=1 

4 

D+X2->DX2 

km 

5 

DX2-)D+X2 

KM 

6 

D+X2->DXS2 

KM 

7 

DXS2->D+X2 

KM 

3 

DX2  +  X2-)DX2X2 

KM 

9 

DX2X2-)DX2+X2 

KM 

10 

DX2+F>DX2+P+4X 

Kf=0.B 

11 

X->A 

Kf:6 

These  are  my  spreadsheet  entries  that  characterize  the  lysis-lysogeny 
pathway  (previously  discussed  in  section  3.7.1). 
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*** **************** OUTPUT  FILE  (lysis-lysogeny  pathway) ****************** 

#  C:\appendix\lysis_lysogeny.odef  Generated  by  JigCell 

#Functions 

#Dependent  species 

DX2X2=(T1  -  DX2  -  DXS2  -  D)/(-(-1.0)) 

P=(T2) 

#Species 

#Independent  Species 

dX2/dt=  -  1*X2  -  1*D*X2  -  1*D*X2  -  1*DX2*X2  +  1*X*X  +  1*DX2  +  1*DXS2  +  1*DX2X2 
dDX2/dt=  -  1*DX2  -  1*DX2*X2  -  ,8*DX2*P  +  1*D*X2  +  1*DX2X2  +  ,8*DX2*P 
dD/dt=  -  1*D*X2  -  1*D*X2  +  1*DX2  +  1*DXS2 
dDXS2/dt=  -  1*DXS2  +  1*D*X2 

dX/dt=  -  2*1*X*X  -  6*X  +  .4  +  2*1*X2  +  4*.8*DX2*P 

dA/dt=6*X 

#Globals 

#Initial  Conditions 

init  X2=. 78934,  0X2=1.23268,  D=l. 56164,  DXS2=1. 23268 

init  X=. 88845,  A=0 

#Constants 

param  T2=1.25,  Tl=5 

#Plot  dependent  species 

aux  DX2X2=DX2X2 

aux  P=P 

done 

*******************0UTPUT  FILE  ( ly s i s_ly sogeny  pathway)  ********************* 


A-3 


Reaction 

Modifiers  and  Constants 

1 

Sig 

A 1  "A3/(A2"A3+A  1  "A3) 

2 

OpSig 

1  -(A1"A3/(A2"A3+A1"A3)) 

3 

kG 

kg  1  *(  1  -  (A  1  /xOmax)) 

4 

->population 

Kh=kG(p  opulation)  *p  opulation 

5 

->luxR 

Kf=T c *(Op Sig(C 0 ,  KcO,  VcO)*Sig(CKP,  KCRP,  VCKP)+b) 

6 

luxR-> 

Kf=l/HRhTA 

7 

luxR-> 

K£=kG(p  opulation) 

8 

-  >luxICD  ABEG 

Kf=Tc*(Sig(C0,  KcO,  VcO)*OpSig(CKP,  KCRP,  VCRP)+b) 

9 

luxICD  ABEG-  > 

K5=1/HBBTA 

10 

luxICD  ABEG-  > 

Kf=kG(p  opulation) 

11 

->LUXR 

Kf=Tl*luxR 

12 

LUXR-> 

Kh=l/Hsp 

13 

LUXR-> 

Kf=kG(p  opulation) 

14 

LUXR+Ai-  >  C  0 

KP=rAIR 

15 

G 

o 

1 

V 

rj 

£ 

Kf=rC0 

16 

->LUXI 

Kf=Tl*luxICDABEG 

17 

LUXI-> 

Kh=l/Hsp 

18 

LUXI-> 

Kf=kG(p  opulation) 

19 

->LUXAB 

KB=H*luxICDABEG 

20 

LUXAB-> 

K^l/Hsp 

21 

LUXAB-> 

Kh=kG(p  opulation) 

22 

->LUXCDE 

Kf=H*luxICDABEG 

23 

LUXCDE-> 

Kh=l/Hup 

24 

LUXCDE-> 

Kf=kG(p  opulation) 

25 

->Ai 

KP=population*(rAII*LIIXI  - 

*** 

rAIR^A^UXR+rC 0 *C 0)+(rAER*LUXR*A  -  rC0*C0) 

26 

Ai-> 

Kh=l/HAI 

27 

C0-> 

Kf=l/Hsp 

28 

C0-> 

Kf=kG(p  opulation) 

Figure  A. 2.  These  are  my  spreadsheet  entries  that  characterize  the  quorum  sensing 

system  that  controls  luminescence  in  V.  fischeri  (previously  discussed 
in  section  3.7.2).  Functions  are  defined  in  lines  1,  2,  and  3  by  corre¬ 
sponding  algebraic  expressions  listed  in  the  Modifiers  mid  Constants 
column. 


A-4 


**********************QUTPUT  FILE  (¥.  fischeri) *************************** 

#  C:\appendix\vfish.odef  Generated  by  JigCell 

#Functions 

Sig(Al,A2,A3)=Al~A3/(A2~A3+Al~A3) 

OpSig(Al , A2 , A3)=l- (A1~A3/ (A2~A3+A1~A3) ) 
kG(Al)=kgl*(l-(Al/xOmax) ) 

#Dependent  species 
#Species 

#Independent  Species 

dluxR_18/dt=  -  l/HRNA*luxR_18  -  kG(popula_26) *luxR_18  + 

(Tc* (OpSig(CO ,  KcO ,  ¥cO)*Sig(CRP,  KCRP,  VCRP)+b)) 
dLUXR_19/dt=  -  l/Hsp*LUXR_19  -  kG (popula_26) *LUXR_19  -  rAIR*LUXR_19*Ai  + 
Tl*luxR_18  +  rC0*C0 

dC0/dt=  -  rC0*C0  -  l/Hsp*CO  -  kG(popula_26) *C0  +  rAIR*LUXR_19*Ai 
dLUXAB/dt=  -  l/Hsp*LUXAB  -  kG(popula_26) *LUXAB  +  Tl*luxICD_25 
dAi/dt=  -  rAIR*LUXR_19*Ai  -  l/HAI*Ai  +  rC0*C0  + 

(popula_26*(rAII*LUXI-rAIR*Ai*LUXR_19+rC0*C0)+ 

(rAIR*LUXR_19*Ai  -  rCO*CO)) 

dLUXCDE/ dt=  -  l/Hup*LUXCDE  -  kG (popula_26) *LUXCDE  +  Tl*luxICD_25 
dLUXI/dt=  -  l/Hsp*LUXI  -  kG(popula_26) *LUXI  +  Tl*luxICD_25 
dluxICD_25/dt=  -  l/HRNA*luxICD_25  -  kG(popula_26) *luxICD_25  + 

(Tc* (Sig(C0 ,  KcO,  VcO)*OpSig(CRP,  KCRP,  ¥CRP)+b)) 
dpopula_26/ dt=kG (popula_26) *popula_26 
#Globals 

#Initial  Conditions 

init  luxR_18=0 ,  LUXR_19=0,  C0=0,  LUXAB=0 
init  Ai=6 ,  LUXCDE=0,  LUXI=0,  luxICD_25=0 
init  popula_26= . 00000015 
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#Constants 


param  Kc0=20,  CRP=10,  VCRP=1,  HRNA=60 

param  Hsp=3600,  rC0=.01,  Tc=0.25,  KCRP=10 

param  Tl=0.25,  b=. 00001,  rAIR= . 000001 ,  Vc0=2 

param  kgl= . 0001925 ,  x0max=.0015,  Hup=600,  rAII=.05 

param  HAI=36000 

#Plot  dependent  species 

done 

**********************0UTPUT  FILE  (¥.  fischeri) ******************************* 
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Reaction 

Modifiers  and  Constants 

->Timer1 

Kf=1 

->Timer2 

Kf=(1/1 2) 

->Alpha 

Kf=0 

-Beta 

KtO 

->CLN2 

Kf=  (ks  n  2 '+ ks  n  2  "*S  B  F)*MAS  S 

CLN2-> 

Kf=kdn2 

->CLB2T 

Kf=(ksb2'+ksb2',*MCM1)*MASS 

CLB2T-> 

Kf=Vdb2 

Vdb2 

kdb2**(HCT1  T-HCT1  )+kdb2”*HCT1  +kdb2'"*CDC20 

->CLB5T 

Kf^(ksb5‘  +ksb5"*MBF)*MASS 

CLB5T-> 

Kf=Vdb5 

Vdb5 

kdb5'+kdb5'*CDC20 

CLN3S 

C  LN  3  MAX*(D  n  3*MAS  S)/ (J  n  3  +  Dn3*MASS) 

->SIC1T 

Kfckscl '  +  kscl  "*SWI5 

SIC1T-> 

Kf=(kd1  cl  +(Vd2c1/(Jd2c1  +SIC1T))) 

->C2 

Kf=kasb2*(CLB2T -  C2)*(SIC1T -(C2+C5)) 

C2-> 

Kf=(kdib2+Vdb2+kd1  cl  +(Vd2c1/(Jd2c1  +SIC1T))) 

->C5 

Kf=kasb5*(CLB5T -  C5)*(SIC1T  -(C2+C5)) 

C5-> 

Kf=(kdib5+Vdb5+kd1  cl  +(Vd2c1/(Jd2c1  +SIC1T))) 

->CDC20T 

Kt(ks20'+ks20'*(CLB2T-C2)) 

CDC20T-> 

Kf=kd20 

->CDC20 

Ktka20*(CDC20T-CDC20) 

CDC20-> 

Kf=(Alpha*Vi20+Beta*((-9.9/1 2)*Timer1  +1 0)+kd20) 

->Vi20 

Kf^O 

->HCT1 

Kf=((kat1  '+kat1  "*CDC20)*(HCT1  T-HCT1  ))/(Jat1  +HCT1 T-HCT1 ) 

HCT1-> 

k1=1;  M1=Vit1;  J1=Jit1 

Vitl 

kitl  '+kit1  ”*(CLN3S+eit1  n2*CLN2+eit1  b5*CLB5+eit1  b2*CLB2) 

Vasbf 

kasbf(CLN2+esbfn3*(CLN3S+BCK2)+esbfb5*CLB5) 

SBF 

G(Vasbf,kisbf+kisbr*CLB2,Jasbf,Jisbf) 

MBF 

G(Vasbf1kisbf+kisbf**CLB21Jasbf,Jisbf) 

MCM1 

G(kamcm*CLB2,kimcm1Jarncm1Jimcm) 

Vd2c1 

kd2cr(ec1  n3*CLN3S+ed  k2*BCK2+CLN2+ec1  b5*CLB5+ec1  b2*CLB2) 

->0RI 

Kf=  ks  o  ri*((C  LB  5T-  C  5) + e  o  ri  b  2*(C  LB  2T  -  C2)) 

0RI-> 

Kf=kdori 

->BUD 

Kfcksbud*(CLN2+CLN3S+ebudb5*(CLB5T  -  C5)) 

BUD-> 

Kf=kdbud 

->SPN 

Kf=ksspn*((CLB2T  -  C2)/(Jspn+(CLB2T-  C2))) 

SPN-> 

Kf=kdspn 

->MASS 

Kfcmu*MASS 

BCK2 

BCK20*MASS 

G 

(2*A1  *A4)/(A2-A1  +A1  *A4+A2*A3+sqrt((A2-A1  +A1  *A4+A2*A3)a2-4*(A2-A1  )*A1  *A4)) 

SWI5 

G(kaswi*CDC20,  kiswi'+kiswi'*CLB2,  Jaswi,  Jiswi) 

CLB2 

CLB2T-C2 

CLB5 

CLB5T-C5 

SIC1 

SIC1T-(C2+C5) 

Figure  A. 3.  These  are  my  spreadsheet  entries  that  characterize  the  budding  yeast 

cell  cycle  (previously  discussed  in  section  3.7.3).  The  single  function 
is  defined  in  line  41.  The  remaining  entries  define  state  variables. 
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>t=H==t==t==t:H==)==t==i«H==t==t==t:H==t=^==i=H==t==i==i==^ OUTPUT  FILE  (budding  yeast)  *************************** 

#  C:\appendix\yeast.odef  Generated  by  JigCell 

#Functions 

G ( A 1 , A2 , A3 , A4) = (2*A1*A4) / (A2-A1+A1*A4+A2*A3+ 

( (A2-A1+A1*A4+A2*A3) ~2-4* (A2-A1) *A1*A4) ~ . 5) 

#Dependent  species 
#Species 

Vdb2=(kdb2 ; )  *  ( (HCT1T) -HCTl)+(kdb2 ’ ’ ) *HCTl+(kdb2J J ,)*CDC20 
aux  Vdb2=Vdb2 

Vdb5= (kdb5 ’ )  +  (kdb5 J ’ ) *CDC20 
aux  Vdb5=Vdb5 

CLN3S= (CLN3MAX) * ( (Dn3) *MASS) / ( ( Jn3)  +  (Dn3)*MASS) 
aux  CLN3S=CLN3S 

Vitl=(kitl ; )+(kitl 5 ’ ) * (CLN3S+(eitln2) *CLN2+(eitlb5) *CLB5+(eitlb2) *CLB2) 
aux  Vitl=Vitl 

Vasbf = (kasbf ) * (CLN2+ (esbf n3) * (CLN3S+BCK2) + (esbf b5) *CLB5) 
aux  Vasbf=Vasbf 

SBF=G(Vasbf , (kisbf 1 )+(kisbf 5 ’ ) *CLB2 , ( Jasbf ) , ( Jisbf ) ) 
aux  SBF=SBF 

MBF=G(Vasbf , (kisbf 1 )+(kisbf 5 ’ ) *CLB2 , (Jasbf) , (Jisbf) ) 
aux  MBF=MBF 

MCMl=G((kamcm)*CLB2, (kimcm) , (Jamcm) , (Jimcm)) 
aux  MCM1=MCM1 

Vd2cl=(kd2cl)*((ecln3)*CLN3S+(eclk2)*BCK2+CLN2+(eclb5)*CLB5+(eclb2)*CLB2) 
aux  Vd2cl=Vd2cl 
BCK2=(BCK20) *MASS 
aux  BCK2=BCK2 

SWI5=G((kaswi)*CDC20,  (kiswi ’ )+ (kiswi ’ 5 ) *CLB2 ,  (Jaswi) ,  (Jiswi)) 
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aux  SWI5=SWI5 


CLB2=CLB2T-C2 
aux  CLB2=CLB2 
CLB5=CLB5T-C5 
aux  CLB5=CLB5 
SIC1=SIC1T- (C2+C5) 
aux  SIC1=SIC1 
#Independent  Species 

dCLB2T/dt=  -  Vdb2*CLB2T  +  ( (ksb2 >+ksb2 ; 5 *MCM1) *MASS) 

dSIClT/dt=  -  ( (kdlcl+(Vd2cl/ (Jd2cl  +  SIC1T) ) ) ) *SIC1T  +  (kscl’  +  kscl,,*SWI5) 
dC5/dt=  -  ( (kdib5+Vdb5+kdlcl+(Vd2cl/ ( Jd2cl+SIC1T) ) ) ) *C5  + 

(kasb5* (CLB5T  -  C5)*(SIC1T  -(C2+C5))) 
dCDC20/dt=  -  ( (Alpha*Vi20+Beta* ( (-9 . 9/ 12) *Timer 1+10) +kd20) ) *CDC20  + 

(ka20* (CDC20T-CDC20) ) 

dHCTl/dt=  (((katl’+katl’ ’ *CDC20) * (HCT1T-HCT1) ) / ( Jatl+HCTIT-HCTI) )  - 
( 1 *HCT 1 * Vit 1 ) / ( J it 1+HCT1 ) 

dBUD/dt=  -  kdbud*BUD  +  (ksbud* (CLN2+CLN3S+ebudb5* (CLB5T  -  C5))) 
dCLB5T/dt=  -  Vdb5*CLB5T  +  ((ksb5’  +ksb5 ’ ’ *MBF) *MASS) 
dCDC20T/dt=  -  kd20*CDC20T  +  ( (ks20 ; +ks20 ; ’ * (CLB2T  -  C2))) 
dORI/dt=  -  kdori*0RI  +  (ksori* ( (CLB5T-C5)+eorib2* (CLB2T  -  C2))) 
dC2/dt=  -  ( (kdib2+Vdb2+kdlcl+(Vd2cl/ ( Jd2cl+SIC1T) ) ) ) *C2  + 

(kasb2* (CLB2T  -  C2)*(SIC1T  -(C2+C5))) 
dSPN/dt=  -  kdspn*SPN  +  (ksspn* ( (CLB2T  -  C2)/(Jspn+(CLB2T  -  C2)))) 
dCLN2/dt=  -  kdn2*CLN2  +  ( (ksn2 ’ +ksn2 5 ; *SBF) *MASS) 
dMASS/ dt=mu*MASS 
d¥i20/dt=0 
dBeta/dt=0 
dTiraerl/dt=l 
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dTimer2/dt=(l/12) 

dAlpha/dt=0 

#Globals 

global  -1  { (CLB2T-C2)  -0.3  }  {MASS  =  f*MASS;  SPN=0.0000;  BUD=0.0000  } 
global  -1  { (CLB2T-C2+CLB5T-C5) -0 . 2  >  {0RI=0.0000  } 
global  1  {0RI-1  }  {Vi20  =10;  Alpha=l  } 

global  1  {SPN-1  }  {Alpha=0 . 0000000000;  Beta=l;  Timerl=0;  Timer2=0  } 
global  1  {Timer2-1  }  {Beta=0 . 0000000000 ;  ¥i20=0.100;  Alpha=l  } 
#Initial  Conditions 

init  CLB2T= . 007 ,  SIC1T=.97,  C5=.03,  CDC20=.0625 
init  HCT1=1,  BUD=.065,  CLB5T=.06,  CDC20T=.2420 
init  0RI=. 015625,  C2=.003,  SPN=.001,  CLN2=.02 
init  MASS=.97,  Vi20=0.100,  Beta=0,  Timer 1=0 
init  Timer2=2,  Alpha=l 
#Constants 

param  Jd2cl=0.05,  HCT1T=1,  kisbf,,=6,  Jisbf=0.01 
param  f=.433,  kd20=0.08,  Jasbf=0.01,  kdlcl=0.01 
param  kisbf)=0.5,  BCK20=.0027,  kaswi=l,  kiswi’=0.3 
param  kiswi,,=0.2,  Jaswi=.l,  Jiswi=.l,  ksspn=0.08 
param  Jspn=0.2,  kdspn=0.06,  mu=0. 005776,  kdori=0.06 
param  ksbud=0 . 3 ,  ebudb5=l ,  kdbud=0 . 06 ,  kamcm=l 
param  kimcm=0.15,  Jamcm=l,  Jimcm=l,  kd2cl=.3 
param  ecln3=20,  eclk2=2,  eclb5=l,  eclb2=0.067 
param  ksori=2,  eorib2=0.4,  kasbf=l,  esbfn3=75 
param  esbfb5=.5,  katl,=.04,  katl,,=2,  Jatl=.05 
param  Jitl=.05,  kitl’=0,  kitl’,=.64,  eitln2=l 
param  eitlb5=.5,  eitlb2=l,  ka20=l,  kdib5=0.05 
param  ks20’ =0.005,  ks20:i  ’=0.06,  kasb2=50,  kdib2=0.05 
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param  kasb5=50,  ksn2’=0,  ksn2,,=.05,  kdn2=0 . 1 

param  ksb2’ =0.002,  ksb2’’=0.05,  kdb2’=0.01,  kdb2,,=2 

param  kdb2,,,=.05,  ksb5’ =0.006,  ksb5,;=0.02,  kdb5’=0.1 

param  kdb5,;=0.25,  CLN3MAX=0 . 02,  Dn3=l,  Jn3=6 

param  kscl’=0.02,  kscl,,=0.1 

#Plot  dependent  species 

done 

OUTPUT  FILE  (budding  yeast)  *************************** 


A-ll 
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